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ABSTRACT 

This  thesis  analyzes  the  dynamic  stability  of  positively  buoyant  submersibles.  Six 
degree -of- freedom  equations  of  motion  are  used  to  compute  steady  state  behavior  with 
motion  restricted  to  the  vertical  plane.  Steady  state  solutions  are  analyzed  for  various 
conditions  of  buoyancy  including  changes  in  (1)  the  amount  of  excess  buoyancy,  (2)  the 
location  of  the  center  of  buoyancy,  (3)  the  location  of  the  center  of  gravity,  as  well  as  (4) 
the  deflection  of  bow  and  stern  planes.  The  equations  of  motion  are  then  linearized  around 
these  steady  state  solutions  to  predict  dynamic  response  in  the  vertical  plane.  The  stability 
of  each  solution  is  then  determined  by  eigen  value  analysis.  The  study  then  expands  the 
analysis  to  include  all  six  degrees  of  freedom  (i.e.,  include  stability  analysis  in  the 
horizontal  plane).  Finally,  numerical  integration  methods  are  used  to  verify  the  results. 
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I.       INTRODUCTION 

A.      PROBLEM   STATEMENT 

Controlling  emergency  ascent  situations  on  submersible  vehicles  such  as  dive  plane 
jam  recovery  is  of  concern  to  the  U.S.  Navy.  In  order  to  control  such  situations,  one  must 
first  be  able  to  predict  the  dynamic  response  of  positively  buoyant  submersibles. 

Dynamic  response  equations  of  motion  descibe  the  manuevering  characteristics  of 
submersible  vehicles  for  six  degrees  of  freedom.  These  equations  assume  constant 
coefficients  for  hydrodynamic  forces  and  moments  approximated  by  zero  frequency  added 
mass  and  damping  terms  plus  the  quadratic  terms  for  drag  forces.  The  constant  coefficients 
vary  for  each  vehicle  and  are  dependent  on  such  things  as  vehicle  body  shape,  location  and 
magnitude  of  vehicle  weight,  location  and  magnitude  of  vehicle  buoyancy,  position  of  bow 
and  stern  planes,  position  of  rudder,  vehicle  speed,  vehicle  mass  characteristics,  vehicle 
hydrodynamic  coefficients,  propeller  rpm  and  control  surface  inputs.  This  thesis  uses  the 
equations  of  motion  and  hydrodynamic  coefficients  for  a  submerged  Mark  IX  Swimmer 
Delivery  Vehicle  (SDV)  developed  by  Smith,  Crane,  and  Summey  [Reference  1:11-16,21- 
31]  to  forecast  the  dynamic  stability  of  a  submersible  in  a  positively  buoyant  condition. 

This  study  begins  by  using  the  six  equations  of  motion  to  compute  the  steady  state 
behavior  of  a  submersible  vehicle  with  motion  restricted  to  the  vertical  plane.  The  steady 
state  solutions  in  the  vertical  plane  are  calculated  for  various  conditions  of  buoyancy 
including  changes  in  the  amount  of  excess  buoyancy,  the  location  ol'  the  center  o( 
buoyancy,  the  location  of  the  center  of  gravity,  as  well  as  the  deflection  of  bow  and  stern 
planes.  The  SDV's  equations  of  motion  are  then  linearized  around  these  steady  state 
solutions  to  predict  dynamic  response  motion  in  the  vertical  plane  for  the  various  conditions 


of  buoyancy.  Several  solutions  are  computed  and  the  stability  of  each  solution  is 
determined  by  eigen  value  analysis.  The  thesis  then  expands  the  analysis  to  include  all  six 
degrees  of  freedom  (i.e.  include  stability  analysis  in  the  horizontal  plane).  Finally, 
numerical  integration  methods  are  used  to  verify  the  results. 

B.      EQUATIONS  OF  MOTION 

The  six  degree  of  freedom  equations  of  motion  for  the  submersible  vehicle  shown  in 
Appendix  A  were  taken  from  Smith,  Crane,  and  Summey  [Reference  1:11-16]. 
Differentiation  with  respect  to  time  is  denoted  by  a  dot  over  the  variable;  e.g.  u  =  ^. 

These  equations  are  referenced  to  a  right-hand  orthogonal  axis  system  fixed  in  the  body 

(vehicle)  as  shown  in  Figure  1.  Since  these  equations  are  in  reference  to  a  body  fixed  axis 

system,  the  Euler  angles  of  pitch  (6),roll  ((])),  and  yaw  (^)  are  used  to  specify  orientation 

with  respect  to  the  inertial  reference  system.  The  rotation  sequence  for  (J),  8  and  ty,  and  the 

Euler  angle  rates  for  <{>,  6  and  H*  shown  in  Appendix  B  were  taken  from  Smith,  Crane,  and 

Summey  [Reference  1:18-20].     Major  variables  and  parameters  as  defined  by    Smith, 

Crane,  and  Summey  [Reference  1:7-10]  are  given  below  : 

1.      Dynamic  Variables 

u,v,w  -  Linear  velocity  components  of  vehicle  with  respect  to  origin  of 

body  axes  system  relative  to  fluid. 

p,q,r  -  Angular  velocity  components  of  vehicle  with  respect  to  body 

axes  system  relative  to  inertial  reference  system. 

X,Y,Z  -  Hydrodynamic  force  components  along  body  axes. 

K,M,N  -  Hydrodynamic  moment  components  along  body  axes. 


2.      Mass  Distribution   Parameters 

m  -  Mass  of  the  flooded  vehicle,  including  the  mass  of  the 

entrained  fluid. 


W  -  Weight  of  the  flooded  vehicle,  including  the  weight  of  the 

entrained  fluid  (  =  g  m  ;  where  g  is  the  acceleration  of  gravity). 


V  -  Displacement  volume  of  the  vehicle. 


B  -  Buoyancy  force  acting  on  the  vehicle  (  =  p  g  V  ).  This  is 

independent  of  the  inertial  mass  distribution  of  the  submersible 
vehicle,  including  whether  or  not  it  is  flooded. 


xG,yG,zG  -  Coordinates  of  the  CG  (center  of  gravity)  in  the  body  axis 

system  (Figure  1).  These  will  depend  on  the  mass 
distribution  of  the  vehicle,  including  the  mass  of  the  entrained 
fluid. 


xB,yB,zB  -  Coordinates  of  the  CB  (center  of  buoyancy)  in  the  body  axis 

system  (Figure  1).  These  are  independent  of  the  mass 
distribution  system,  but  may  vary  with  the  addition  or  removal 
of  external  appendages. 


I  ,1  ,1  -  Moments  of  inertia  about  the  body  system  axes,  including  the 

entrained  fluid. 


I    ,1    ,1  -  Products  of  inertia  about  the  body  system  axes,  including  the 

xy    xz    yz  J    J  '  6 

entrained  fluid. 


3.      Remaining  Parameters 

p  -  Mass  density  of  fluid  medium 


Reference  length  used  to  nondimensionalize  the  hydrodynamic 
coefficients. 


b(x),  h(x)  -  Width  and  height  of  vehicle  in  its  xy  and  xz  planes, 

respectively,  at  location  x  measured  in  the  body  axes  system 
(Figure  1).  These  quantities  are  required  in  the  integrals 
defining  the  crossflow  forces  and  moments  in  the  equations  of 
motion,  and  are  tabulated  within  the  Steady  State  Computer 
Program  (Appendix  C). 


x       ,x    .,  -  Sternplane,  bowplane  and  rudder  deflection  angles  in  radians 

(Figure  1). 
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II.     SYSTEM  SOLUTIONS  IN  THE  VERTICAL  PLANE 

A.      GENERAL 

In  the  steady  state  condition,  the  submersible  will  have  reached  constant  linear  and 
angular  velocities.  Therefore,  the  body  fixed  linear  accelerations  (li  ,  v,  w)  and  the  body 
fixed  angular  accelerations  (p  ,  q  ,  r)  will  be  zero.  Similarly,  the  vehicle  will  have  reached 
a  constant  angle  of  pitch  (6)  making  its  derivatives  (6)  equal  to  zero.  Since  this  analysis  is 
restricted  to  steady  state  solutions  in  the  vertical  plane,  the  angle  of  roll  (<t>)  and  its 
derivative  ($)  will  be  zero  (in  Chapter  IV,  the  case  where  the  angle  of  roll  (<{>)  is  180 
degrees  will  be  discussed  during  the  numerical  integration  analysis).  The  angle  of  yaw  (^) 
and  its  derivative  (^)  will  likewise  be  zero  due  to  the  vertical  plane  restriction.  It  should  be 
noted  that  had  this  analysis  not  been  restricted  to  the  vertical  plane,  steady  state  yaw  Qr) 
would  not  necessarily  be  zero,  thereby  allowing  the  angular  velocities  (p,  q,  r)  to  be  non- 
zero. However,  since  this  analysis  was  restricted  to  the  cases  where  (j>,  6  ,  ^'  ,  p,  q  and  r 
are  all  zero,  the  equations  of  motion  for  six  degrees  of  freedom  for  the  steady  state 
condition  can  be  reduced  to: 

•  Longitudinal  (Surge)  Equation  of  Motion: 

( W  -  B)  sin  6  =  [  Xvv  v2  +  Xww  w2  +  X^  uv6r  +  uw  (  Xw6s  6S  +  Xw6b  6b)] 

X6s6s  6s    +  X6b6b  6b    +  X6r6r  6r  I  +  u2xprop 

•  Lateral  (Sway)  Equation  of  Motion: 

-  (W-B)  cos  6  sincf>=  [  Yv  uv  +  Yvw  vw  +  Y^r  u2  6r 

/"xnose 

[  CDy  h(x)  (v)2  +CDz  b(x)  (w)2  ]  ^L-  dx 

'"tail 


Normal  (Heave)  Equation  of  Motion 

-  (W-B)  cos  0  cos  <j>  =  Z^  uw  +  ZvV  v2  +  u2  (zfts  6S+Z&b  65) 


/*xnose 
■'xtail 


:CDyh(x)(v)2+CDzb(x)(w)2]u^x)dx 


•  Roll  Equation  of  Motion: 

-  (yG  W  -  yg  B)  cos  6  cos  (j)  +  (zqW  -  zgB)  cos  0  sin  ty 
=  Kv  uv  +  Kvw  vw]  +  u2  Kpr0p 

•  Pitch  Equation  of  Motion: 

(xq  W  -  xg  B)  cos  6  cos  (p  +  (zq  W  -  zg  B)  sin  6  = 

[  Mw  uw  +  Mvv  v2  +  u2  (M6s  6s  +  M6b  6b)] 

Jrxnose 
[  CDy  h(x)  (v)2  +CDz  b(x)  (w)2  ]  JJ_  x  dx 
xtail  Cl 

•  Yaw  Equation  of  Motion: 

-  (xq  W  -  xg  B)  cos  6  sin  $  +  (yQ  W  -  yg  B)  sin  8  = 

Nv  uv  +  Nvw  vw  +  N6r  u2  6r^  +  u2  Nprop 

xnose 

[  CDy  h(x)  (v)2  +  CDz  b(x)  (w)2  ]  ^-  x  dx 


/"xno 
•^tail 


B.      CONDITIONS 

1 .      Defining  Additional  Terms 
a.      Excess  Buoyancy,   bB 

Excess  buoyancy  is  defined  as  6B  =  B  -  W  where  B  is  the  submersible's 
total  buoyancy  and  W  is  the  submersible's  total  weight. 


b.  Longitudunal  Center  of  Buoyancy,  xqB 

The  longitudinal  center  of  buoyancy  is  defined  as  xqd  =  xg  "  xB  wnere 
x<3  is  the  longitudinal  center  of  gravity  with  respect  to  the  body  fixed  axis  and  Xg  is  the 

longitudinal  center  of  buoyancy  with  respect  to  the  body  fixed  axis. 

c .  Vertical  Center  of  Buoyancy,     zG  B 

The  vertical  center  of  buoyancy  is  defined  as  Zqt>  =  Zq  -  Zt>  where  Zq  is 
the  longitudinal  center  of  gravity  with  respect  to  the  body  fixed  axis  and  Zg  is  the 
longitudinal  center  of  buoyancy  with  respect  to  the  body  fixed  axis  (zqd  is  assumed  to  be 

positive). 

2 .      Assumed   Conditions 

a.  Lateral  Centers  of  Gravity,  yQ  ,  and    Buoyancy,  y™ 

The  lateral  center  of  gravity  and  the  center  of  buoyancy  are  assumed  to  be 
on  the  same  centerline  plane  (yG  =  yB  =  0). 

b .  Propeller  Speed,n     (revolutions  per  minute) 

The  propeller  speed  is  assumed  to  be  zero  (n  =  0). 

c.  Propeller   Coefficients,   Kprop  and  Nprop 

From  Smith,  Crane,  and  Summey  [Reference  1:30],  the  propeller 
coefficients  are  zero  (Kprop  =  Nprop  =  0). 

C.     REVISING  THE  EQUATIONS  OF  MOTION 

Using  the  term  for  vertical  center  of  buoyancy  (zq-d),  the  expression  (zqW  -  ZgB) 

may  be  written  as  (zqrW  -  Zr>5B).  Similarly,  using  the  term  for  longitudinal  center  of 
buoyancy  (xq^)  the  expression  (xqW  -  XgB)  may  be  written  as  (*QgW  -  Xg8B).  Also, 
the  term  u2Xpr0p  may  be  written  as: 

u2Xprop  =  u2CDo(ri2-l)=u2CDo[(^^)2-l]  =  CdoAV  -  Cdou2 


where  A  is  a  constant. 

Since   the  shall  speed   (n)  is  zero,   the  expression  may   be   further   reduced   to 
u2Xprop  =  -  Cdqu2.    Substituting  these  expressions  plus  the  term  for  excess  buoyancy 

(6B)  and  the  assumed  conditions  revises  the  equations  of  motion  for  the  six  degree  of 
freedom  system  as  follows: 

•  Longitudinal  (Surge)  Equation  of  Motion: 

-6B  sin  0  =  [  Xvv  v2  +  Xww  w2  +  Xy6r  uv6r  +  uw  (  Xw6s  6S  +  Xw6b  6b) 

+  lu2  (  X6s6s  6s2  +  X6b6b  6b2  +  X6r6r  63i  "  C™  u" 

•  Lateral  (Sway)  Equation  of  Motion: 
6B  cos  0  sincf)=    Yv  uv  +  Yvw  vw  +  Y^   u2  6r 


I 


xnose 

v       dx 


CDyh(x)v2+CDzb(x)w2JIJ^ 
xtail  C 


Normal  (Heave)  Equation  of  Motion: 

6B  cos  6  cos  <p  =  !  Zw  uw  +  Zvv  v2  +  u2  (z^s  os+Z^b  65) 


I 


xnose 


ICDyh(x)v2+CDzb(x)w2!u^x) 
xtail  C 


dx 


Roll  Equation  of  Motion: 
(zQgW  -  zg6B)  cos  6  sin  <j)=  :  Kv  uv  +  Kvw  vwj 

Pitch  Equation  of  Motion: 

(XGB  W  "  XB  °B)  cos  8  cos  (f)  +  (zqb  W  -  zg  6B)  sin  6 
Mw  uw  +  Mvv  v2  +  u2  (m6s  6s  +  M6b  6^). 


fxnose 
+  [  CDy  h(x)  v2  +CDz  b(x)  w2  1  u^  x  dx 

•'xtail 


•  Yaw  Equation  of  Motion: 

(-XGB  W  +  xg  6B)  cos  0  sin  ty  =  [  Nv  uv  +  Nvw  vw  +  N^r  u2  6r 

/"xnose 

[  CDy  h(x)  v2  +  CDz  b(x)  w2  ]  _L  x  dx 

-'xtail 

These  six  equations  only  have  five  unknowns:  u,v,w,6,  and  <J>.  Therefore,  two  of 
these  equations  must  be  dependent  and  additional  conditions  are  required  in  order  to  make 
the  number  of  equations  equal  the  number  of  unknowns. 

D.     ADDITIONAL  CONDITIONS 

The  next  condition  to  be  applied  to  the  vehicle  is  that  the  rudder  will  remain 

centerlined,  that  is  6r  =  0.    Since  the  solutions  of  interest  are  those  in  which  the  vehicle 

remains  within  the  vertical  plane,  it  can  be   further  specified  that  the  linear  velocity  in  the 

transverse  direction  (v)  equals  zero.  Recalling  that  the  angle  of  roll  (<j))  has  been  previously 

assumed  to  equal  zero,  the  trigonometric  functions  of  <j>  can  be  identified  as  sin<{>  equals  zero 

and  cos<|>  equals  unity.    Substituting  these  quantities  back  into  the  equations  of  motion, 

three  of  the  six  equations  of  motion  (sway,  roll,  and  yaw)  yield  trivial  solutions.     In 
addition,  the  cross-flow  velocity  term  (U  f)  for  the  heave  and  pitch  equations  can  be 

reduced  to: 

Ucf(x)  =  [(v  +xr)2  +  (w  -  xq)2j°-5  =  [w2]0'5  =  |w| 
since  v,  r,  and  q  are  zero.    Furthermore,  since  CD  is  constant,  it  can  be  taken  outside  the 

integral.  Therefore,  the  three  remaining  equations  can  be  wriiten  as: 

•  Longitudinal  (Surge)  Equation  of  Motion: 

AP««fl  Xwww2+uw(Xw6s6s  +  Xw6b6b) 

od  sin  o  =  |  /  j  2\ 

+  u2lX6s6s6s    +  X6b6b6b  )-CD0U2_ 
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Normal  (Heave)  Equation  of  Motion: 

6B  cos  B  =  Z^  uw  +  u2  (z6s  os+Z6fe  6b)  -  CDz  I 


xnose 


b(x)    ,wy    dx 

w 


x,  the  three 


•  Pitch  Equation  of  Motion: 

(XGB  w  "  XB  ^B)  cos  8  +  (zqb  W  -  zg  6B)  sin  6  = 

pnose 

Mw  uw  +  u2  (M6s  6S  +  M6b  6b)+  CDz  b(x)  f7  x  dx 

i  w 

'xlail 

E.  VERTICAL  PLANE  EQUATIONS  OF  MOTION 

By  defining  the  terms  Aw  as  the  b(x)  dx  and  xA  as  -£-         b(x)  x  dx, 

entail  •'xiai! 

remaining  equations  defining  motion  in  the  vertical  plane  can  be  written  as: 

•  Longitudinal  (Surge)  Equation  of  Motion: 

-  6B  sin  6  =  Xww  w2+  uw  (  Xw6s  6S  +  Xw6b  6b) 

+  u2  (  X6s6s  6s2  +  X6b6b  6b2)  -  Cdo  u2 

•  Normal  (Heave)  Equation  of  Motion: 

6B  cos  6  =  Z^  uw  +  u2  (z6s  6S+Z6b  6b)  -  cDz  w  lwl  Av 

•  Pitch  Equation  of  Motion: 

(XGB  W  -  xg  6B)  cos  6  +  [zqb  W  -  zg  6B]  sin  6  = 
Mw  uw  +  u2  (m6s  6s  +  Mb[)  6b)+  CDz  w  |w|  xAA,v 

F.  COMPUTER  PROGRAM  DEVELOPMENT 

Taking  these  three  equations  which  describe  motion  in  the  vertical  plane,  solving  the 
first  two  for  sine  6  and  cosine  6,  respectively;  and  dividing  all  three  through  by  u-  yields: 
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Longitudinal  (Surge)  Equation  of  Motion: 


sin  6  =      1 
u2  6B 


Xww(^)2+(^)(Xw6s6s  +  Xw6b6b 
6s6s  6s    +  X6b6b  6b  )  -  Cdo 


Normal  (Heave)  Equation  of  Motion: 


cos 


-  =  ™  I  *»  u  +  (Z6s  6s+Z6b  6b)  "  CDz  ™  A. 
6B  u^ 


•  Pitch  Equation  of  Motion: 

(xGBw-xB6B)e^e  +  (ZGBW.ZB6B)^iiLe  = 

uz  uz 

M\v  u  +  (M6s  6s  +  M6b  6b)+  CDz  ^M  xAAw 

Now,  defining  the  quantity  ~  as  w',   and  substituting  w'   back  into  the  three 
equations: 

if  w  is  positive: 

•  Longitudinal  (Surge)  Equation  of  Motion: 

sin  6  =      1    |  Xww  (w'l)2+  (w,)(  Xw6s  6s  +  Xw6b  6b 
u2  6B 


(  X6s6s  6s    +  X6b6b  6bl  -  Cdo 


•  Normal  (Heave)  Equation  of  Motion: 

^^[Zw^'  +  ^s  &s+Z6b  6b)  -CDz(w')2A 
u"        oB 

•  Pitch  Equation  of  Motion: 

(xGBW-xB&B)£QS-e   +  (zgbw-zB6B)^9  = 
u2  uz 

Mw  (Wl)  +  (M6s  6s  +  M6b  6b)  +  CDz  (w1)2  xAK 
however,  if  w  is  negative: 

•  Longitudinal  (Surge)  Equation  of  Motion: 


■V. 


sin  6  =      1 
u2  6B 


XWw(w,i)2+(w')(Xw6s6s  +  Xw6b6b) 
+  (  X6s6s  6s2  +  X6b6b  6b^) "  Cdo 
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Normal  (Heave)  Equation  of  Motion: 

co^6  =  J    :z^(w.)  +  (z^  6S+Z6b  6b)  +CDz(wfAw] 
u2        5B 

Piteh  Equation  of  Motion: 

(xgbW-xB6B)^^  +  (zGBw-zB^B)>:il1^  = 
u2  u2 

Mw  (W)  +  (M6s  6s  +  M5b  6b)  -  CDz  (w')2  xAAv 


sin  8  cos  Q_ 

Substituting  the  equations  ior       2      an(^       2      int0  me  pitch  equation  yields  the 

following  expressions: 


if  w  is  positive: 

(xGBW-xB6B)  -^-[^(w'J  +  fz^  6S+Z6b  6b)  -CDz(wf  A,v 
6B 
(         W         ,R)    1       Xww(wf+(w',)K6s6s  +  Xw6b6b 
-  UgBw-zBoB)  /  2  2^ 

6B 


X6s6s  6s"  +  X&b&b  5b  J  -  CDo 


Mw  (w')  +  (M6    6s  +  M6b  6b)  +  CDz  (w')2  xAAvv 


and  if  w  is  negative: 

(xGBW-xB6B)  J    ^(w'j  +  fz^  6S+Z6b  6b)  +CDz(w')2Aw 


6B 


(Xww(w')2+(w',)(Xw6s6s  +  Xw6b6b) 
-  (zGB  W  -  zB  6B)  -i-  /  2      v  *  2\    „ 

^Bi       +  \  x6s6s  6s    +  X6b6b  6b  ) "  CD0 

MwK)  +  (M6s  6S+M6b  6b)-CDz  (wfxA 


Rearranging  these  two  equations  to  get  them  into  the  form:   A(w')-  +  B(w')  +  C  =  0,  the 
expressions  become: 
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Lww 


(wf 


if  w   is  positive: 

Cj)z  xAAw  +  (zqb  w  -  ZB  6B)  —  x 

6B 

+  (xGB  w  "  *B  6B)  ~  CDz  K 
6B 

Mw  +  (zGB  W  -  zB  6B)  ■*-  (  X^  6S  +  Xw6b  6b) 

6B 
-(xGbW-xB6B)  -1-Zw 
6B 


(W) 


(zgb  W  -  zB  6B)  M  X6s6s  6s2  +  X6b6b  6b2  -  CD0) 


6B 


+  Ks  6s  +  M6b  6b)-(xGBW-xB6B)  -Uz6s  6S+Z6b  6b) 

oB 


=  0 


J 


and  if  w  is  negative: 

-  Cdz  xaAw  +  (zgb  W  -  zb  6b)  -      Xww 

6B  (w')2 

-  (xGB  W  -  xB  oB)  J-  CDz  Aw 

6B  J 

Mw  +  (zGB  W  -  zB  SB)  -L  (  Xw6s  6S  +  Xw6b  6b) 

6B 


(w'l) 


-(xcBW-xEdBl-^Zw 
6B 

(zgb  W  -  zB  6B)  -U  X6s6s  6s2  +  X6b6b  6b2  -  CD0) 
6B 
+  (M6s  6S  +  M6„  6b)-(xGBW-xB6B)  -Uz6s  6S+Z6„  6y) 

oB 

These  quadratic  expressions  were  then  solved  using  the  equation: 


=  0 


w"  = 


-B±  \B2-4AC 

2A 


where: 


B  =6B  Mw  +  (zGB  W  -  zB  6B>  (  X^  6S  +  Xw6b  6b) 

-(*GB  W-xb&b)  2V 
C  =  (zgb  W  -  zB  6B)  (  X6s6s  6S2  +  X6b6b  6b2  -  CD0) 

+  6B(M6s  6S  +  M6b  6b)-(xGBW-xB6B)   (z6s  6S+Z6b  6b) 

if  w  is  positive: 

A  =  6B(CDz  xAAw)  +  (zgb  W  -  zb  6B)     Xww  +  (xgb  w  -  *B  6B)   CDz  A*.. 
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and  if  w  is  negative: 

A  =  -  oB(CDz  xA\j)  +  (zGB  W  -  zB  6B)     Xww  -  (xGB  W  -  xB  6B)    CDz  Au 

The  value  of  8  was  determined  using  the  computed  values  of  w'  and  the  equation  for 

tangent  6: 

sin  8 
tan  8  =  -"?- 


cos  0 
u2 


However,  the  value  of      2     varied  depending  on  the  value  of  w',  which  lead  to  two 
possible  solutions: 

Equation  for  tan  8   if  w'   is  Positive: 

tan  8  =    -^w^-Hl^sV  Xw6b  6b)  -(  X6s6s  6S2  +  X6b6b  6b2) kC^ 

Zw(w')++  (Z6s  6S+Z6b  6b)  -CDz  A^w'KI 


Equation  for  tan  8   if  w'   is  Negative: 

tan  8  =    -Xww  'w'»2  -  'w''  (  *w6s  5s  +  *w6b  5b)  -  (  Xfe6s  6S2  +  X6bftb  6b2)  +  CD0 

ZVJ(W')++  (Z6s  6S+Z6b  5b)  +CDz  ^(w'Jw'l 

In  either  case,  the  value  of  u2  was  computed  using  the  expression  derived  from  the 
surge  equation  of  motion: 


u2 6B  sin  8 

-  Xww  (w'i)2  -  (w1)  (  Xw6s  6S  +  Xw6b  6b)  -  (  X6s6s  6S2  +  X&b6b  6b2)  +  CD0 

This  leads  to  two  possible  solutions  for  u  (i.e.  u  =  ±  \  u2).  The  value  of  w  was 
computed  using  w  =  u  (w').  Combining  the  two  possible  solutions  of  u  with  the  two 
possible  values  for  w'  derived  from  the  quadratic  expression  lead  to  four  possible 
combinations  of  solutions  for  u  and  w.  The  computer  program  which  calculates  these  four 
possible  solutions  is  contained  in  Appendix  C.  It  is  an  interactive  program  designed  to 
allow  the  operator  to  select  the  amount  of  excess  buoyancy  as  a  percentage  of  vehicle 
weight  (6B),  the  deflection  of  dive  planes  in  degrees  (os),  the  ratio  of  bow  planes  to  dive 
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weight  (6B),  the  deflection  of  dive  planes  in  degrees  (6S),  the  ratio  of  bow  planes  to  dive 
planes  (6b/6s),  the  location  of  xGB  and  xB  from  body  fixed  axis  origin  as  a  percentage  of 
length,  and  the  location  of  zGB  and  zB  from  body  fixed  axis  origin  in  feet. 

G.     STEADY  STATE  RESULTS 

Figures  2,  3,  and  4  show  typical  steady  state  solutions  for  surge  velocity,  heave 
velocity,  and  pitch  angle  as  a  function  of  dive  plane  angle.  The  two  cases  shown  vary  the 
location  of  the  longitudinal  center  of  buoyancy;  for  case  (a):   xGB  =  -  1  %  of  the  vehicle 

length  (L),  and  for  case  (b):  xGB  =  +  1  %  L.  The  following  parameters  were  the  same  for 

both  cases:  excess  buoyancy  ,  6B  =  2  %  of  the  vehicle  weight  (W);  deflection  of  bow 
planes,  6b  =  0;  location  of  horizontal  and  vertical  centers  of  buoyancy,  xR  =  zR  =  0;  and 
location  of  vertical  center  of  gravity,  zGB  =  0.1  feet. 

All  runs  developed  four  solutions.  For  two  of  the  solutions,  the  magnitude  of  the 
surge  velocitities  were  large  while  the  magnitudes  of  the  associated  heave  velocities  were 
relatively  small.  This  has  been  descibed  by  Booth  [Ref  2:  297]  as  "predominantly  forward 
motion".  The  other  two  solutions  had  small  surge  velocity  magnitudes  and  large  heave 
velocity  magnitudes.  Booth  [Ref  3:  346]  referred  to  this  type  of  motion  as  "nearly  vertical 
ascents".  The  positive  or  negative  nature  of  the  velocities  are  associated  with  the  value  of 
pitch  angle.  Positive  heave  velocities  are  associated  with  pitch  angles  greater  than  90 
degrees.  That  is  the  submersible  would  be  ascending  in  a  belly  up  orientation.  Although 
this  steady  state  analysis  computes  four  possible  solutions,  it  gives  no  indication  as  to 
which  of  the  solutions  are  stable  (if  any).  Dynamic  response  and  stability  criteria  will  be 
discussed  in  the  next  chapter. 
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Figure  2.     Steady   State   Vertical   Plane  Solutions  for  Surge   Velocity   (u) 
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Figure  3.     Steady  State  Vertical  Plane  Solutions  for  Heave  Velocity  (w) 
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Figure  4.     Steady  State  Vertical  Plane  Solutions  for  Pitch  Angle  (6) 
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III.  DYNAMIC  STABILITY 


A.  GENERAL 

The  first  portion  of  this  chapter  uses  the  six  degree  of  freedom  equations  of  motion 
along  with  the  Euler  angle  rate  equations  for  the  derivatives  of  the  angles  of  pitch  and  roll 
(6  and  <|))  to  predict  the  dynamic  stability  when  movement  is  restricted  to  the  vertical  plane. 
These  equations  of  motion  are  then  linearized  around  the  vertical  plane  steady  state  nominal 
points  computed  in  chapter  II.  Eigen  value  analysis  is  then  used  to  compute  the  stability  of 
each  solution.  The  second  part  of  the  chapter  expands  this  analysis  to  include  all  six 
degrees  of  freedom  and  uses  the  same  steady  state  nominal  points  to  predict  the  dynamic 
stability  when  motion  is  not  restricted  to  the  vertical  plane. 

B.  RESTRICTING  MOTION  TO  THE  VERTICAL  PLANE 

Since  motion  is  restricted  to  the  vertical  plane,  the  body-fixed  transverse  veloctiy  (v) 
and  its  derivative  (v)  are  zero.  The  angles  of  roll  (<p)  and  yaw  (^  ),  and  their  derivatives  (<|> 
and  ^)  are  also  zero.  We  will  continue  to  assume  that  the  lateral  center  of  gravity  and  the 
lateral  center  of  buoyancy  are  on  the  same  centerline  plane  (yQ  =  ye  =  0)>  and  the  rudder  is 
centerlined  (or  =0). 

C.  LINEARIZED  VERTICAL  PLANE  EQUATIONS  OF  MOTION 

Substituting  these  values  into  the  original  equations  of  motion  (Appendix  A),  yields 
trivial  results  for  three  of  the  six  equations:  Lateral  (sway),  Roll,  and  Yaw.  The  remaining 
equations  reduce  to  the  following  form: 
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Longitudinal  (Surge)  Equation  of  Motion: 
|  m  -  Xu:  li  +  [m  zG  }  q  =  [  -CD0  +  X6s6sos2  +X6b6b6b2  j  u2 
+  [Xw6s6s  +  Xw6b6bJ  uw  +  Xq6s6s  +  Xq6b6bj  uq+  [Xw]  w- 
+  [Xwq  -  m]  wq+  'Xqq  +  mxG  q2  -  (W  -  B)  sin  9 

Normal  (Heave)  Equation  of  Motion: 

m  -  ZvJ  w  -  m  xq  +  Zq'  q  =  j  ZAs6s  +  Z6bftb  u2 
+  Zw]  uw  +  i  mj  uq+  [mzc1  q2 


- 

J  SI: 


[ CDz  b(x)(w-xq)2]  ^-f-  dx  +  (W  -  B)  cos  6 

ucfW 


Pitch  Equation  of  Motion: 

[mzo]  u  - 1  Mw  +  m  xG]  w  +  [  Iy  -  Mq  q  =  M^&s  +  M^bbj  u2+  [Mw]  uw 


+  i Mq  -  mxo  uq-  mzc  wq 

-  (xGW  -  xBB)  cos  6  -  (zgW  -  zBB)  sin  8 


Jf  xdosc 
x  tail 


x  dx 


These  three  equations  which  describe  motion  restricted  to  the  vertical  plane  are  functions  of 
four  variables  (u,w,q,0)  plus  their  deriatives  (u,w,q,6).  The  equations  of  motion  and  the 
Euler  angle  rate  equation  for  8  were  linearized  using  the  following  generalized  procedure: 


Bll(i,l)u  +  Bll(i,2)w  +  Bll(i,3)q  +  Bll(i,4)8  = 


rafi 


du 


u  + 


u, 


dfil 


dw 


\dt 


w  + 

wo     1.  c*q 


q  + 


dfi 


*       ^0J0O 


t) 


where  the  Bll(i)'s  are  the  constants  associated  with  the  derivatives  of  the  variables,  the 
functions  fj  represents  the  right  hand  side  of  the  nonlinear  equations,  and    i  =   1  to  4 
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identifies  the  three  equations  of  motion  (surge,  heave,  and  yaw)  plus  the  Euler  angle  rate 
equation  for  6.  The  partial  derivatives  were  computed  as  follows: 

•         Partial  Derivatives  of  the  Longitudinal  (Surge)  EOM: 

-±  =  2(-CDo  +  X6s6s6s2  +X6b6b6b2)u0  +  (xw6s6s  +  Xw6b6b)wo  =  Al  1(1,1) 


^  =  2Xww0  +  (Xw6s6s  +  Xw6b6b)u0  =  Al  1(1,2) 
ow 


rff  /  \ 

y1  =  (Xwq  -  m)wo+  (Xq6s6s  +  Xq6b6b)u0  =  Al  1(1,3) 


^  =(W-B)  cos  6o  =  Al  1(1,4) 
•        Partial  Derivatives  of  the  Normal  (Heave)  EOM: 

f2   =  Z..W0  +  2  (Z6s6s  +  Z6b6b)u0  =  Al  1(2,1) 

OU 


~ 2  =  uoZw  -  2CDzA>0|  =  All(2,2) 
dw 


—*■  =  u0(Zq  +  m)  +  20)^x^01  =  Al  1(2,3) 
oq 


^  =  -(W-B)sine0  =  All(2,4) 

ae 


Note  on  differentiation  procedure:  The  cross-flow  velocity  term  (Ucf)  was 
reduced  to  Ucf(x)  =  [(v  +xr)2  +  (w  -  xq)2]05  =  [w2]05  =  |w  -  q| .  Allowing 
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the  integral  term  in  the  heave  equation  to  be  defined  as  L,  where: 
h  =  Q)7  (         b(x)(w-xq)|w  -  xq|  dx 


■ 

/xiail 


d(w  -  q) 


«xnose  m\ 

/ 1         b(x)(w-xq)sign  (w  -  xq)  dx  =  2CD/  I 

Atail  Jxti 


=  2Cdz  I         b(x)(w-xq)sign  (w  -  xq)  dx  =  2CD/  J         b(x)  |w  -  xq|  dx 


dl.      I       d\ 


dw       \  d(w 


*xnosc 

V")  ( ^^) =  2CDidWo1 1    b(x)  dx = 2Cdz  |wo1  Aw 


J*  x  nose 
xb(x)dx  =  -2CDz|w0|xAAv 
xiail 


•        Partial  Derivatives  of  the  Pitch  EOM: 


^  =  2  (M5S6s  +  M^bjuo  +  Mw  w0  =  Al  1(3,1) 

du 


^  =  Mvv  u0  +2u0CDzAwX>0|  =  Al  1(3,2) 


di-x 

-^  =  (Mq  -  mxC})u()  -(mzG)w(|  -  2CDz  I^w0|  =  Al  1(3,3) 
dq 


^  =  (xGW  -  xBB)  sinBo  -  (zGW  -  zBB)  cos  60  =  Al  1(3,4) 
d9 


Note  on  differentiation  procedure:  Again  the  cross-flow  velocity  term  (U,f)  was 
reduced  to  |w  -  q|.  Allowing  the  integral  term  in  the  pitch  equation  to  be  defined 
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as  L,  where: 

Mxaosc 

1 2  =  Cqz  I         b(x)(w-xq)|w  -  xq|  x  dx 


J  XV, 


d(w  -  q) 


J  XV, 


-  2Cdz  I         b(x)  |w  -  xq|  x  dx 


Jaxnose 
b(x)xdx  =  2CD;|w0|xAA. 
x  tail 


=(^V)(^)-2C^'  ^h(x)dx  =  -2CDz,w0|IA 


dq      \  d(w  -  q)   /  \      dq 

Partial  Derivatives  of  the  Euler  Angle  Rate  Equation  for  6: 

^  =  0  =  Al  1(4,1) 

du 


^i=0  =  Al  1(4,2) 
dw 


^i=  1  =  Al  1(4,3) 
dq 


^  =  0  =  Al  1(4,4) 

ae 


The  constants  corresponding  to  the  derivatives  of  the  four  variables  (i.e.  the  left  hand 
side  of  the  four  equations)  arc  as  follows: 
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Constants  Associated   with  Derivatives  for  the  Longitudinal 
(Surge)   EOM: 


u    =>    m-XL1  =  Bll(U) 
q    =*    mzG  =  Bll(l,3) 

Bl  1(1,2)  =  Bl  1(1,4)  =  0 

•  Constants  Associated  with  Derivatives  for  the  Normal 

(Heave)   EOM: 

w    =>  m-Zw  =  Bl  1(2,2) 

q    ^-(mxG  +  Zq)  =   Bl  1(2,3) 

Bl  1(2,1)=  Bl  1(2.4)  =  0 

•  Constants  Associated  with  Derivatives  for  the  Pitch  EOM: 

u    =>  mzCj  =  Bl  1(3,1) 

vv    =>  -(Mw  +  mxG)=Bl  1(3,2) 


q    =*>    Iy-M^  =  B1 1(3,3) 
B  11(3,4)  =  0 

•       Constants  Associated  with  Derivatives  for  the  Euler  Angle  Rate 
Equation  for  0  : 

9    =>   1  =B1 1(4,4) 

Bl  1(4,1)  =  Bll(4,2)  =  Bl  1(4,3)  =0 
These  expressions  can  be  arranged  in  a  matrix  format  to  form  the  linearized  equations 
of  motion  in  the  vertical  plane  about  the  nominal  steady  state  points.  The  matrix  format  is 
as  follows: 


Bll  xXl  =  All  x  XI 


where: 


All  = 


Bll  = 


Al  1(1,1) 

All(l,2) 

Al  1(1,3) 

All(l,4) 

All(2,l) 

Al  1(2,2) 

Al  1(2,3) 

Al  1(2,4) 

Al  1(3,1) 

Al  1(3,2) 

Al  1(3,3) 

Al  1(3,4) 

Al  1(4,1) 

Al  1(4,2) 

Al  1(4,3) 

Al  1(4,4) 

Bl  1(1,1) 

Bll(l,2) 

Bll(l,3) 

Bll(l,4) 

Bll(2,l) 

Bl  1(2,2) 

B  11(2,3) 

B  11(2,4) 

Bll(3,l) 

Bl  1(3,2) 

Bl  1(3,3) 

Bl  1(3,4) 

Bll(4,l) 

Bl  1(4,2) 

Bl  1(4,3) 

Bl  1(4,4) 
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and 

u 

v  1  w 

XI  =  i 

q 
_  e 

D.  VERTICAL  PLANE  COMPUTER  PROGRAM  DEVELOPMENT 

The  matrix  format  of  the  linearized  equations  of  motion  was  used  in  the  computer 
program  shown  in  Appendix  D  to  predict  the  dynamic  stability  of  the  vehicle  with  motion 
restricted  to  the  vertical  plane.  The  program  is  interactive  in  that  it  allows  the  operator  to 
select  which  of  the  four  data  files  from  the  steady  state  analysis  (Chapter  II)  will  be  used  to 
define  the  nominal  points  for  the  linearization  process.  An  eigen  system  subroutine  was 
used  to  find  eigen  values  and  eigen  vectors.  The  program's  output  was  called  the  degree  of 
stability  and  only  shows  the  largest  real  part  of  all  eigen  values.  The  stability  criteria  is 
such  that  the  degree  of  stability  must  be  negative  in  order  for  the  solution  to  be  stable. 

E.  VERTICAL  PLANE  DYNAMIC  RESULTS 

Of  the  four  possible  steady  state  solutions  computed  in  Chapter  II,  only  one  solution 

from  each  case  yielded  stable  characteristics.  There  were  some  cases  in  which  none  of  the 

solutions  were  stable  for  certain  ranges  of  parameters.   The  general  trend  of  the  linearized 

dynamic  results  are  fairly  well  demonstrated  by  the  two  cases  discussed  in  Chapter  II. 
Recalling  the  parameters  of  these  cases:  6B  =  2%,  ratio  of  bow  to  stern  planes  (o./6s)  =  0, 

ZGB  =  ^-1  ^eet'  ar,d  xr  =  ZB  =  ^'  ^ne  f'rst  case  P'acec*  the  longitudinal  center  of  gravity  aft 
of  the  longitudinal  center  of  buoyancy  (xGB  =  -  1  %),  and  the  second  case  placed  the 
longitudinal  center  of  gravity  forward  (xGB  =  +  1  %).  Once  again  dive  plane  deflection 
angle  (6  )  was  varied  from  -  20  to  +  20  degrees  for  both  cases.  Figure  5  shows 
longitudinal  velocity  (u)  as  a  function  of  dive  plane  angle  (6  ).    Case  One  (xGB  =  -  1  %) 
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showed  predominantly  forward  motion,  while  case  two  (xGB  =  -  1  %)  yielded  nearly 

vertical  ascents.  Figure  6  shows  vertical  motion  (w)  as  as  a  function  of  dive  plane  angle 
(&s).   The  results  concur  with  Figure  5,  case  one  shows  very  little  vertical  motion  while 

case  two  demonstrates  a  larger  value.  It  is  interesting  to  note  that  vertical  motion  (w)  for 
case  one  is  positive  for  dive  plane  angles  between  -  20  and  -4  degrees.  The  case  one  values 
of  6  shown  in  figure  7  for  dive  plane  angles  between  -  20  and  -4  degrees  concur  with  this 
observation.  The  values  of  0  greater  than  90  degrees  indicate  the  submersible  is  ascending 
in  the  belly  up  position.  The  stability  in  the  vertical  plane  is  shown  in  figure  8.  Degree  of 
stability  (e)  is  shown  as  a  function  of  dive  plane  angle  (6s).   This  figure  shows  both  cases 

to  be  stable  in  the  vertical  plane. 

12 


10- 


Figure  5.     Stable  Vertical   Plane  Solutions  for  Surge  Velocity 
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-0.5 


Figure  6.     Stable  Vertical   Plane  Solutions  for  Heave  Velocity 
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Figure  7.     Stable  Vertical  Plane  Solutions  for  Pitch  Angle  (6) 
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Figure  8.     Degree  of  Stability  in  the  Vertical  Plane 

F.      LINEARIZATION  OF  FULL  EQUATIONS  OF  MOTION 

Referring  back  to  the  complete  equations  of  motions  shown  in  Appendix  A,  these  six 
equations  which  describe  motion  for  the  six  degree  of  freedom  system  are  functions  of 
eight  variables  (u,v.w,p.q.r,<p.G)  plus  their  deriatives  (u,v,w,p,q,r,<j>,0).  The  six 
equations  of  motion  (Appendix  A)  along  with  the  Euler  angle  rates  (Appendix  B)  were 
linearized  as  folows: 

bjlu  +  bj2v  +  bj3w  +  bj4p  +  bj5q  +  bj6r  +  bj7<p  +  bj89  = 


dgj  ;  d&  '  dg; : 

-s-    u+  -5L    v+  -5L     w  + 
dujuo       Ldv|Vo       [dwJWo 


^gj 


dgj 


dr 


ro 


r+  -SL    4>+  -M\    e 

to  j<j)0     L  ^e  Je0 


L  dP  jPo       L  dq  Jqo 

where  the  bj's  are  the  constants  associated  with  the  derivatives  of  the  variables,  the 
functions  g;  represents  the  right  hand  side  of  the  nonlinear  equations,  and    j  =  1  to  8 
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identifies  the  six  equations  of  motion  plus  the  Euler  angle  rate  equations  for  <J>  and  6. 
respectively.  The  partial  derivatives  of  these  equations  were  computed  as  follows: 

•         Partial  Derivatives  of  Longitudinal   (Surge)   EOM: 

-f1  =(Xw6s  6s  +  Xw6b  6b)w„  +2(X6s6s  6s~  +  X5h6b  6b  )u0  -  2CD„u0  =  all 
du 


-^  =  2Xww0  +  (Xw6s  6s  +  Xw6b  6b)u0  =  a  13 
dw 


-p  =  -mw0  +  Xwqw0  +  (Xq6s  6s  +  Xqftb  6b)u0  =  a  15 


^-  =  -(W-B)cose  =  al8 

d6 


*gi=0=al2  ^  =  0  =  al4  *8Uo  =  al6  ™  =  0  =  a17 

dv  dp  dr  d(f> 


Partial  Derivatives  of  Lateral  (Sway)  EOM: 


;-:  =Yvu0  +  Yvww0  -  CDAv  kol  =  a22 
dv 


p  =mw0  +  Yp  u()  +Yvvpw0  =  a24 
dp  '  ' 


~^-2  =  -  muo  +Yr  u()+  Ywr  w0  -  CD/Awx^w0|=  a26 
dr 


dg2  =  (W  -  B)  cos6  =  a27 
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dg2 
du 


=  0  =  a21 


^g2-n- 


dw 


=  0  =  a23 


dg2 


=  0  =  a25 


^  =  0=a28 

de 


Note  on  differentiation  procedure:  The  cross-flow  velocity  term  (U  f)  is  given 
by:  Ucf(x)  =  [(v  +xr)2  +  (w  -  xq:)2j  '  .  The  integral  term  in  the  sway  equation 


was  defined  as  L,  where: 


u  = 


xnose 


r 


CDy  h(x)  (v+xr)2  +CDz  b(x)  (w-xq j2  j  *-±XL  dx 


mose 


'^))dx 


dv 


dv 


(y-uf)  +  I^(^)  =  c-A»w«2 


|w0| 


Ucf  -  (v  +  xr)  —  Ucf 
dv 


Ui 


=  CDzw02  ]  -z  Aw  =  CDz  |w0|  Av 
w02 


ai3    /  ai 


dr 


du 


dr 


( ^l +I  f ( Vr1) =  Cdz AwW°2  "V  x  Urf " (v  +  xr)^F 

\    Ucf  7      dr  \    Ucf  /  ]j2f  [  dr 


=  CDz  Aww02  Xa|V!0'   =  Cdz  Awxa  |w0| 
w02 


ik_  =  ( JL\  | v  +  xr)|  +1 1_  I v  +  xr))  =  c 

dw        \  dw    /  \     Ucf    /       dw  \     Ucf 


dh  =  I  J]_\  I  v_+xO  |  +I  A_  I  v  +  xr)j  =  0 
dq       \dq   }\     Ucf    )      dq\     Ucf    / 
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•         Partial  Derivatives  of  Normal  (Heave)  EOM: 


-^  =  Zw  w0  +  2  (Z6s  6s  +  Zfth  6b')u0  =  a31 
du 


-^-  =  Z,v  u()  -  2  CDz  A,v  |w0|  =  a33 
dw 


-f-  =  m  u0  +  Zq  u0  +  2  CDz  Aw  x^w0|  =  a35 
dq  ^ 


?&  =.  (W  -  B)  sin  e0  =  a38 
d6 


^=0=a32  ^3=0=a34  ^  =  0  =  a36  ™  =  °  =  a37 

dv  dp  dr  d(f) 


Note  on  differentiation  procedure:  The  integral  term  in  the  heave  equation 
was  defined  as  L,  where: 

/"xnose 
I4  =  CDy  h(x)  (v+xr)2  +CDz  bfo)  (w-xq)2  ]  ^^  dx 

Jxtail  '  « 


dv       \3v  l\      Urf     I       3v\      Urf  I       (Ucfp 

/^      *         ->    1  dUcf      „  .  dUcf      „ 

=  -    Cdz  Awwq2-^4t  wq~  =  0  because    — ^    =  0 

w0" 


ci 


dv  dv 
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dl4      /   dl  W  w-xq\  d    /  w  -  xq 

aw       \  dw  )  \     Ucf    /        dw  \     Ucf 

=  2CDzA^w-xq)^a  +  I^9 

ucf       ucf2  L 

=  2  CD/  Aww02-L  +cDz  A^wo2-1 

|W(,I 


Ucf-(w  -xq) 


aUef 
dw 


w02. 


|w0|  -  w0 


Wq 

]wol 


=  2CDzAJw0| 


because    ^  =  |ucf(x)  =  [(v  +xr)2  +  (w  -  xqfj0  52(w  -  xq)  =  g 


w 


dl4  _  /   dl  \  /  w  -  xq 


dq        \  dq 


ra-i-ra-^ 


x  Ucf  -  (w  -  xq}— 


ducf 


=  CDz  Aww02^^-  =  CDz  Awxa|w()| 
w02 


dr 


di4 


=  0 


Partial  Derivatives  of  Roll  EOM: 


^  =  Kv  uo  +  Kvvv  w0  =  a42 
dv 


d24 

-f5-  =  -  m  zG  w0  +  Kp  uo  +  KWp  w0  =  a44 

dp  F       l 


d2i 

~-  -    m  zq  uo  +  Kr  uo  +  Kwr  w0  =  a46 
dr 

^-4  =  -  (  zG  W  -  zB  B)  cos90  =  a47 

d(j> 

^-  =  0  =  a41            ^4  =  0=a43 

dg4 

=  0  = 

a45 

du 


dw 


dq 


^  =  0  =  a48 

de 
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Partial  Derivatives  of  Pitch  EOM: 


^  =  Mw  w0  +  2  (NV  6s  +  H,b  6b)u0  =  a51 
du 


-r^5  =  Mw  u0  +  2  CDz  AwX^wol  =  a53 
dw 


-f-  =-  mxG  u0  -  mzG  w0  +  H>  u0  -  2  CDz  I^w()|  =  a55 
dq  ^ 


^  =  (xGW  -  xBB)  sin  0O  -  (zGW  -  zBB)  cos  80  =  a58 
d6 


^  =  0  =  a52  ^=0  =  354  ^>=0  =  a56  ~  =  °  =  a57 

dv  dp  dr  d(p 


Note  on  differentiation  procedure:  The  integral  term  in  the  pitch  equation 
was  defined  as  1^,  where: 


I5  = 

tail 


/"Xnosc 

[  CDy  h(x)  (v+xr)2  +CDz  b(x)  (w-xq)2  ]  ^-  x  dx 

•'xtail  Cf^ 

rn0Sem/w-xq;\       .       fXn0SCT         , 

i   (i)yixdx=xri  i4xdx 


:tail  /xtail 

dl5        dU 


dw         dw 


dh       dl4 


dq         dq 


xA=  2CDzAwxA|w0| 


xA  =  CDz  Awxa2|w0|  =  CDz  IaIwqI 


dis       dU  „  dl.s       dl4  n 

-7-^-  =  t-^  xa  =  0  - —  =  ^ —  xA=  0 

dv        dv  dr        dr 
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Partial  Derivatives  of  Yaw  EOM: 

dg6 


dv 


=  Nv  uo  +  Nvw  w0  -  CDz  Aw  xA  |w0|  =  a62 


-%£-=  m  xG  w0  +  Np  uo  +  Nwp  w0  =  a64 

dp  p  v 


-^---  m  xG  u0  +  Nr  u0  +  Nwr  w0  -  CDz  IA  |w0|  =  a66 
dr 


^=    (xG  W  -  xB  B)  cos0o  =a67 
d(j) 


*&  =  0  =  a61  *6>  =  0  =  a63  ^  =  0=a65  ^  =  0  =  a68 

du  dw  dq  d8 


Note  on  differentiation  procedure:  The  integral  term  in  the  yaw  equation 
was  defined  as  L,  where: 

Jf  xnose 
[  CDy  h(x)  (v+xr)2  +CDz  b(x)  (w-xq)2  ]  ^r  x  dx 
xtail  C 

Jrxnose  /*xnose 

x.,  (I«)xH  ,  l3Xdx 
xtail  •'xtail 

^6  Cll3  /-I        I         I    A 

t —  =  t —  xA  =  CDz  w0  Aw  XA 
dv        dv 


T^"  =  ™  XA  =  CDz  A^xa2  |w0|  =  CDz  IA|w0| 
dr        dr 


al3  xA=o  i^=i^XA  =  o 


dw         dw  dq       dq 
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Partial  Derivatives  of  Euler  Angle  Rate  for  (p 


^7  =  1  =  a74 

dp 

^  =  tan  60  = 

ar 

a76 

^7=0=a71 

du 

av 

=  0  =  a72 

c*g7 

aw 

=  0  =  a73 

^7  =  0  =  a75 
dq 

dg? 

=  B=a77 

dg7 

ae 

=  0  =  a78 

Partial  Derivatives  of  Euler  Angle  Rate  for  6  : 


Hfi*  =  1  =  a86 

aq 

^-=0=a81 

au 

^  =  0  =  a82 

av 

^  =  0  =  a83 

aw 

^  =  0  =  a85 

aq 

^gs  _  n  _  ,o. 

^  =  0  =  a87 

^  =  0  =  a88 

dr  a<p  ae 


The  constants  corresponding  to  the  derivatives  of  the  eight  variables  (i.e.  the  left  hand 
side  of  the  eight  equations)  are  as  follows: 
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•  Constants  Associated   with  Derivatives  for  the  Longitudinal 

(Surge)   EOM: 

u      =>     (m-Xu)=bll  q      =>     (mzG)  =  bl5 

bl2  =  bl3  =  bl4  =  bl6  =  bl7  =  bl8  =0 

•  Constants  Associated  with  Derivatives  for  the  Lateral  (Sway)  EOM: 

v     =>    (ra-Yv)  =  b22  p     =>    (-mzo  -  Yp)  =  b24 

r      =>     (mxo  -  Yf)  =  b26 
b21  =  b23  =b25  =  b27  =  b28  =  0 

•  Constants  Associated  with  Derivatives  for  the  Normal  (Heave) 

EOM: 

w      =>     (m  -  Z^,)  =  b33  q      =>     (-  mxG  -  Z^)  =  b35 

b31  =  b32  =  b34  =  b36  =  b37  =  b38 

•  Constants  Associated  with  Derivatives  for  the  Roll  EOM: 

v      =>     (-mzG  -  Kv)  =  b42  r      =>     (  -  Kr)  =  b46 

p      =>     (Ix-Kp)  =  b44 
b41  =  b43  =  b45  =  b47  =  b48 

Constants  Associated  with  Derivatives  for  the  Pitch  EOM: 

u      =>     (mzo)=b51  w      =>     (- mxG  -  Mv,)  =  b53 

q      =*     (Iy-Mq)  =  b55 
b52  =  b54  =  b56  =  b57  =  b58 

•  Constants  Associated  with  Derivatives  for  the  Yaw  EOM: 

v      =>     (mxG-Nv)  =  b62  p      =>     (-Np)=b64 

f   =>  (Iz  -  Nr)  =  b66 
b61  =  b63  =  b65  =  b67  =  b68  =  0 
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•  Constants  Associated   with   Derivatives  for  the  Euler  Angle  Rate 

Equation  for  (p  : 

4      =>     1  =  b77 

b71  =  b72  =  b73  =  b74  =  b75  -  b76  =  b78  =  0 

•  Constants  Associated  with  the  Derivatives  for  the  Euler  Angle 

Rate  Equation  for  6  : 

e     =>    i  =  b88 

b81  =  b82  =b83  =  b83  =  b84  =  b85  =  b86  =  b87  =  0 

These  expressions  can  be  arranged  in  a  matrix  format  to  define  the  dynamic  equations 

of  motion  for  the  six  degree  of  freedom  system  linearized  about  the  nominal  steady  state 

points.  The  matrix  format  is  B  x  X  =A  x  X,  where: 

T  all     al2    al3    al4    al5    al6    al7    al8 

a21  a22  a23  a24  a25  a26  a27  a28 

a31  a32  a33  a34  a35  a36  a37  a38 

a41  a42  a43  a44  a45  a46  a47  a48 

a51  a52  a53  a54  a55  a56  a57  a58 

a61  a62  a63  a64  a65  a66  a67  a68 

a71  a72  a73  a74  a75  a76  a77  a78 

.  a81  a82  a83  a84  a85  a86  a87  a88 


A    = 


all 

0 

al3 

0 

al5 

0 

0 

a  18 

0 

a22 

t) 

a24 

0 

a26 

a27 

0 

a31 

0 

a33 

0 

a35 

0 

0 

a38 

0 

a42 

0 

a44 

0 

a46 

a47 

aO 

a51 

t) 

a53 

0 

a55 

0 

0 

a58 

0 

a62 

0 

a64 

0 

a66 

a67 

0 

0 

0 

0 

1 

0 

a76 

(i 

0 

0 

0 

0 

0 

1 

0 

(i 

0 
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B  = 


bll  bl2  bl3  bl4  bl5  bl6  517  bl8 

b21  b22  b23  b24  b25  b26  b27  b28 

b31  b32  b33  b34  b35  b36  b37  b38 

b41  b42  b43  b44  b45  b46  b47  b48 

b51  b52  b53  b54  b55  b56  b57  b58 

b61  b62  b63  b64  b65  b66  b67  b68 

b71  b72  b73  b74  b75  b76  b77  b78 

b81  b82  b83  b84  b85  b86  b87  b88 


and 


bll 

0 

0 

0 

bl5 

0 

0 

0 

0 

b22 

0 

b24 

0 

b26 

0 

0 

0 

0 

b33 

0 

b35 

0 

0 

0 

0 

b42 

0 

b44 

0 

b46 

0 

0 

b51 

0 

b53 

0 

b55 

0 

0 

0 

0 

b62 

0 

b64 

0 

b66 

0 

0 

0 

0 

0 

0 

0 

0 

1 

0 

0 

0 

0 

0 

0 

0 

0 

1 

x= 


xl 

u 

x2 

V 

x3 

w 

x4 

p 

x5 

q 

x6 

r 

x7 

♦ 

„x8^ 

e 
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These  matrices  may  be  reordered  such  that  they  will  be  of  the  form: 


B 1 1      0 
0     B22 


XI 

.  X'2  - 


All      0 

0     A22 


[XI 
X2 


This  is  accomplished  by  rewriting  the  X  matrix  such  that  XI  is  the  same  matix  used 
in  the  vertical  plane  analysis: 


X  = 


xl 

u 

x3 

w 

x5 

q 

x8 

0 

x4 

p 

x7 

♦ 

x2 

V 

x6 

r 

XI 
X2 


The  A  matrix  is  restructured  into  a  matrix  containing  four  4x4  matrices  with 
theA12  and  A21  matrices  containing  only  the  zero  element. 


A  = 


all 

al3 

al5 

al8 

0 

0 

0 

0 

a31 

a33 

a35 

a38 

0 

0 

0 

0 

a51 

a53 

a55 

a58 

0 

0 

0 

0 

0 

0 

1 

0 

0 

0 

l) 

0 

All   0 

0 

0 

0 

0 

a44 

a47 

a42 

a46 

0  A22  - 

0 

0 

0 

0 

1 

0 

0 

a  76 

0 

0 

0 

0 

a24 

a  27 

a22 

a26 

0 

0 

0 

0 

a64 

a67 

a62 

a66 
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Similarly,  B  is  restructured  into  a  matrix  containing  four  4x4  matrices  with 
the  B12  and  B2 1  matrices  containing  only  the  zero  element. 


B  = 


)11 

0 

bl5 

0 

0 

0 

0 

0 

0 

b33 

b35 

0 

0 

0 

0 

0 

)51 

b53 

b55 

0 

0 

0 

0 

0 

0 

0 

0 

1 

0 

0 

0 

0 

_iBll   0 

0 

0 

0 

0 

b44 

0 

b42 

b46 

L  0  B22 

0 

0 

0 

0 

0 

1 

0 

0 

0 

0 

0 

0 

b24 

0 

b22 

b26 

0 

0 

0 

0 

b64 

0 

b62 

b66  _ 

As  discussed,  the  XI  matrix  established  to  descibe  the  linearized  dynamic  stabilty  in  the 
vertical  plane  is  the  same  as  the  XI  matrix  within  X.  In  addition,  the  All  and  Bll 
matrices  from  the  vertical  plane  analysis  are  also  identical  to  those  in  X  .  That  is: 


All  = 


Al  1(1,1)  All(l,2)  All(l,3)  Al  1(1,4) 

Al  1(2,1)  Al  1(2,2)  Al  1(2,3)  Al  1(2,4) 

Al  1(3,1)  Al  1(3,2)  Al  1(3,3)  Al  1(3,4) 

Al  1(4,1)  Al  1(4,2)  Al  1(4,3)  Al  1(4,4) 


a44  a47  a42  a46 

1  0  0  a76 

a24  a27  a22  a26 

a64  a67  a62  a66 


Bll  = 


Bll(l,l)  Bll(l,2)  Bll(l,3)  Bll(l,4) 

Bl  1(2,1)  B1J(2,2)  Bl  1(2,3)  Bll(2,4) 

Bll(3,l)  Bl  1(3,2)  Bl  1(3,3)  Bl  1(3,4) 

Bl  1(4,1)  Bll(4,2)  Bl  1(4,3)  B  11(4,4) 

The  A22,  B22,  and  X2  matrices  represent  the  additional  equations  necessary  to  describe 

the  linearized  motion  for  all  six  degrees  of  freedom;  henceforth  referred  to  as  the  horizontal 


b44  0   b42  b46 

0  10  0 

b24  0   b22  b26 

b64  0   b62  b66 
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plane  contributions.  The  eigen  function  for  the  six  degree  of  freedom  model  is  computed 
by  taking  the  determinant  as  follows: 


det 


A11-XB11  0 

0  A22-XB22 


=  0  =>  (det!  All  -XBll])(det[A22-XB22])  =  0 


Since  the  eigen  function  will  be  the  product  of  these  two  determinants,  the 
resulting  eigen  values  will  merely  be  the  union  of  the  vertical  plane  eigen  values  and  the 
horizontal  plane  eigen  values.  The  significance  of  reducing  the  eigen  value  calculation  from 
an  8  by  8  matrix  problem  to  two  4  by  4  matrix  problems  is  not  in  the  computation  time 
saved.  But  rather  in  fact  that  now  the  horizontal  and  vertical  stabilities  have  been  separated 
and  identified. 

G.     COMPUTER  PROGRAM  DEVELOPMENT 

The  matrix  format  of  the  linearized  dynamic  response  equations  associated  with 
motion  in  the  horizontal  plane  was  added  to  the  computer  program  developed  previously 
(Appendix  D).  Once  again,  the  program  is  interactive  in  that  it  allows  the  operator  to  select 
which  of  the  four  data  files  from  the  steady  state  analysis  (Chapter  II)  will  be  used  to 
define  the  nominal  points  for  the  linearization  process.  An  eigen  system  subroutine  was 
used  to  find  eigen  values  and  eigen  vectors.  Two  outputs  were  added  to  the  program. 
First,  the  degree  of  stability  in  the  horizontal  plane,  and  next  the  degree  of  stability  of  both 
planes  (i.e.,  the  union  of  the  vertical  and  horizontal  degrees  of  stability).  Reminder: 
degree  of  stability  must  be  negative  in  order  for  the  solution  to  be  stable. 

H.     DYNAMIC  STABILITY  SOLUTIONS 

Continuing  with  the  same  cases  from  part  E,  the  horizontal  stability  for  the  two  cases 
are  shown  in  figure  9.    Case  two  (*CR  =  +  1  %)  is  stable  in  the  horizontal  plane  for  all 
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values  of  dive  plane  angle  (6S).  Whereas,  case  one  (xGB  =  -  1  %)  is  unstable  in  the 
horizontal  plane  for  dive  plane  angle  (6S)  between  -20  and  -  9  degrees.  From  figure  7,  this 
corresponds  to  values  of  6  greater  than  140  degrees.  This  indicates  that  the  vehicle  is 
stable  (even  in  the  horizontal  plane)  for  values  of  0  greater  than  90  degrees.  A  submersible 
with  an  angle  of  pitch  greater  than  90  degrees  will  have  a  negative  metacentric  height  and 
will  therefore  be  statically  unstable.  However,  the  results  shown  in  figures  7  and  9  indicate 
that  the  vehicle  will  remain  dynamically  stable  is  such  a  condition.  This  'inverted 
pendulum'  type  stability  will  be  further  investigated  during  the  numerical  integration 
analysis  (Chapter  IV)  to  see  if  hydrodynamic  and  drag  forces  on  the  vehicle  can  actually 
cause  this  to  occur.  Figure  10  shows  the  combined  stabilities  for  the  horizontal  and  vertical 
planes.  It  should  be  noted  that  in  general  the  horizontal  plane  dictated  stabilty  for  the  cases 
considered. 
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Figure  9.     Degree  of  Stability  in  the  Horizontal  Plane 
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Figure  10.     Degree  of  Stability  in  Both  (Horizontal  and  Vertical)  Planes 
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Because  the  real  part  of  the  computed  eigen  values  must  be  negative  in  order  for  a 
solution  to  be  stable,  it  is  easy  to  identify  unstable  solutions.  However,  the  measure  of 
stability  for  a  stable  solution  is  harder  to  quantify.  The  magnitudes  shown  thus  far  for 
degree  of  stability  seem  very  small.  Does  this  mean  that  the  solutions  are  not  very  stable? 
In  an  attempt  to  answer  this  question,  the  values  for  degree  of  stability  for  a  neutrally 
buoyant  submersible  (6B  =  0  or  W  =  B)  were  computed.  Figure  11  shows  the  degree  of 
stability  in  (a)  the  horizontal  plane  and  (b)  the  vertical  plane  as  a  function  of  surge  velocity. 
This  figure  shows  that  the  magnitude  of  the  values  for  degree  of  stability  for  this  neutrally 
buoyant  case  are  indeed  of  the  same  order  of  magnitude  as  the  positively  buoyant  cases 
discussed  earlier. 
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Figure  11.     Degree  of  Stability  as  a  Function  of  Surge  Velocity  (u)  for  a 

Neutrally   Buoyant  Submersible 
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IV.    SIMULATIONS  AND  DISCUSSION   OF  RESULTS 

A.     DYNAMIC  STABILITY   ANALYSIS   RESULTS 

The  dynamic  stability  analysis  included  in  this  Chapter  considers  stability  for  both 

planes  (i.e.  the  analysis  no  longer  breaks  stability  down  by  horizontal  or  vertical  plane). 

As  mentioned  in  Chapter  III,  horizontal  plane  stability  generally  dictates  the  stability  of  the 

vehicle. 

1 .      Variations  in  Longitudinal  Center  of  Gravity  (xGB) 

Figures  12  through  15  show  how  changing  the  longitudinal  center  of  gravity 
(xGB)  effects  the  dynamic  response  of  the  submersible.    For  these  cases,  the  amount  of 

excess  buoyancy  (6B)  is  two  percent  of  wieght  (W),  the  bow  plane  deflection  angle  (6.)  is 

zero,  the  vertical  center  of  gravity  (zGB)  is  0.1  feet,  and  the  longitudinal  and  vertical  centers 

of  buoyancy  (xB  and  z^)  are  zero.  The  longitudinal  center  of  gravity  (xCR)  is  varied  from 

-  2  to  +  2  percent  of  vehicle  length.  When  the  longitudinal  center  of  gravity  (xrR)  is  greater 

than  zero,  there  are  stable  solutions  for  the  full  range  of  dive  plane  angles.   However,  this 
is  not  true  when  the  longitudinal  center  of  gravity  (xGB)  is  less  than  zero.    The  range  of 

stable  solutions  is  restricted  when  xGB=  -  0.5  and  -  1.0. 
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(a)  xGB   =  0,  0.5,  1.0,   1.5  and  2.0   % 
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(b)   xGB  =  0,  -  0.5,  -  1.0,  -  1.5  and  -  2.0  % 


Figure  12.     Stable  Surge  Velocity  (u)  Solutions  for  Variations  in  xGB 
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(b)   xGB   =  0,  -  0.5,  -  1.0,  -  1.5  and   -  2.0   % 


Figure   13.     Stable   Heave  Velocity   (w)   Solutions  for  Variations  in   xGB 
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(b)  xGB  =  0,  -  0.5,  -  1.0,  -  1.5  and  -  2.0  % 


Figure  14.     Stable  Angle  of  Pitch  (8)   Solutions  for  Variations  in  xGB 
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Figure   15.     Degree  of  Stablity  for  Variations  in  xGB 
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2.      Variations  in  the  Amount  of  Excess  Buoyancy  (8B) 

Figures  16  through  19  show  how  changing  the  amount  of  excess  buoyancy 

(8B)  effects  the  dynamic  response  of  the  submersible.    For  these  cases,  the  bow  plane 
deflection  angle  (8.)  is  zero,  the  longitudinal  center  of  gravity  (xGB)  is  -  0.5  percent  of 

vehicle  length,  the  vertical  center  of  gravity  (zGB)  is  0.1  feet,  and  the  longitudinal  and 

vertical  centers  of  buoyancy  (xfi  and  zg)  are  zero.  The  amount  of  excess  buoyancy  (8B)  is 

varied  from  0.5  to  2.9  percent  of  weight  (W).    The  only  case  where  there  are  stable 
solutions  for  the  full  range  of  dive  plane  angles  is  when  buoyancy  (8B)  is  0.5. 
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(b)  6B   =  2.5,   2.6,  2.7,  2.8,  2.9   % 


Figure   16.      Stable  Surge   Velocity   (u)   Solutions  for  Variations  in   6B 
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Figure   17.     Stable  Heave  Velocity  (w)   Solutions  for  Variations  in  6B 
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Figure  18.     Stable  Angle  of  Pitch  (6)   Solutions  for  Variations  in  6B 
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(b)  6B   =  2.5,  2.6,  2.7,  2.8,  2.9    % 


Figure  19.     Degree  of  Stablitj   for  Variations  in  6B 


56 


3.      Variations   in    Vertical   Center  of  Gravity   (zGR) 

Figures  20  through  23  show  how  changing  the  vertical  center  of  gravity  (zrR) 

effects  the  dynamic  response  of  the  submersible.  For  these  cases,  the  amount  of  excess 
buoyancy  (6B)  is  two  percent  of  weight  (W),  the  bow  plane  deflection  angle  (6.)  is  zero, 
the  longitudinal  center  of  gravity  (xGB)  is  -  0.5  percent  of  vehicle  length,  and  the 
longitudinal  and  vertical  centers  of  buoyancy  (xB  and  zR)  are  zero.  The  vertical  center  of 
gravity  (zGB)  is  varied  from  0.05  to  0.30  feet.  The  general  trend  shown  in  these  figures  is 
that  the  larger  values  for  vertical  center  of  gravity  (zGB)  have  stable  solutions  over  a  larger 
range  of  dive  plane  angles. 
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Figure  20.     Stable  Surge  Velocity  (u)  Solutions  for  Variations  in  z 
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Figure  21.     Stable  Heave  Velocity   (w)  Solutions  for  Variations  in  zGB 
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Figure  22.     Stable  Angle  of  Pitch  (6)  Solutions  for  Variations  in  zGB 
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Figure  23.     Degree  of  Stability  for  Variations  in  z 
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4.      Variations  in  the  Longitudinal  Center  of  Buoyancy   (xB) 

Figures  24  through  27  show  how  changing  the  longitudinal  center  of  buoyancy 
(xB)  effects  the  dynamic  response  of  the  submersible.  For  these  cases,  the  amount  of 
excess  buoyancy  (oB)  is  two  percent  of  weight  (W).  the  bow  plane  deflection  angle  (6.)  is 
zero,  the  longitudinal  center  of  gravity  (xGB)  is  -  0.5  percent  of  vehicle  length,  the  vertical 
center  of  gravity  (zGB)  is  0.1  feet,  and  the  vertical  center  of  buoyancy  (zB)  is  zero.  The 
longitudinal  center  of  buoyancy  (Xg)  is  varied  from  -  9  to  +  2  percent  of  vehicle  leneth. 

The  general  trend  shown  in  these  figures  is  that  the  positive  values  for  longitudinal  center 
of  buoyancy  (Xg)  lend  to  have  stable  solutions  for  positive  dive  plane  angles.  On  the  other 
hand,  the  negative  values  for  longitudinal  center  of  buoyancy  (x&)  tend  to  have  stable 
solutions  for  negative  dive  plane  angles. 
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Figure  24.     Stable  Surge  Velocity  (u)  Solutions  for  Variations  in  xB 
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Figure  25.      Stable   Heave   Velocity   (w)   Solutions  for  Variations  in   xfi 
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Figure  26.     Stable  Angle  of  Pitch  (6)   Solutions  for  Variations  in  x, 
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Figure  27.     Degree  of  Stability  for  Variations  in  x. 
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5.      Variations  in  Bow  Planes  Deflection  Angle  (6b) 

Figures  28  through  31  show  how  a  non-zero  bow  planes  deflection  angle  (6.) 

effects  the  dynamic  response  of  the  submersible.  For  these  cases,  the  amount  of  excess 
buoyancy  (6B)  is  two  percent  of  wieght  (W),  the  longitudinal  center  of  gravity  (xGB)  is 
-  0.5  percent  of  vehicle  length,  the  vertical  center  of  gravity  (zGB)  is  0.1  feet,  and  the 
longitudinal  and  vertical  centers  of  buoyancy  (xB  and  Zg)  are  zero.    The  bow  planes  are 

given  a  deflection  of  -  20  degrees.  The  significance  of  the  results  shown  in  these  figures  is 
that  for  certain  dive  plane  angles  (6s  =  -3  to  -12  degrees)  there  are  two  stable  solutions. 
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Figure   28.      Stable   Surge   Velocity   (u)   Solutions   for  a   Non-zero  Bow   Plane 

Angle  (6b  =  -  20  degrees) 
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Figure  29.     Stable   Heave   Velocity   (w)   Solutions   for  a   Non-zero  Bow   Plane 

Angle  (6.    =  -  20  degrees) 
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Figure  30.     Stable  Pitch  Angle  (6)  Solutions  for  a  Non-zero  Bow  Plane 

Angle  (&b  =  -  20  degrees) 
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Figure  31.     Degree  of  Stability   for  a   Non-zero  Bow  Plane  Angle 

(6b  =  -  20  degrees) 

B.      SIMULATIONS   USING   NUMERICAL  INTEGRATION  METHODS 

The  linearized  dynamic  response  results  were  verified  by  simulations  using  numerical 
integration  of  the  full  six  degrees  of  freedom  equations  of  motion  for  the  swimmer  delivery 
vehicle  (SDV).  Figure  32  shows  a  plot  of  angle  of  pitch  (6)  versus  time  for  the  center  of 
gravity  forward  of  the  center  of  buoyancy  case  (xGB  =  +  1)  with  a  dive  plane  angle  (&s)  of 

-  15  degrees.  The  dotted  line  shows  the  linearized  results  from  figure  7,  the  solid  line 
shows  the  numerical  integration  results.  The  steady  state  results  of  the  numerical 
integration  method  match  the  linearized  dynamic  results  exactly. 
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Figure  32.     Numerical  Integration  Solution   for  Angle  of  Pitch   (8)   when 
Center  of  Gravity  is  Forward  (xGB  =  +  1%) 

Figure  33  shows  a  plot  of  angle  of  pitch  (9)  versus  time  for  the  center  of  gravity  aft 
of  the  center  of  buoyancy  case  (xGB  =  -  1)  with  a  dive  plane  angle  (o  )  of  -  15  degrees. 

Again  the  dotted  line  shows  the  linearized  results  from  figure  7.  the  solid  line  shows  the 
numerical  integration  results.  And  once  again,  the  steady  state  results  of  the  numerical 
integration  method  match  the  linearized  dynamic  results  exactly.  However,  this  linearized 
dynamic  result  was  for  the  vertical  plane  only,  the  horizontal  plane  stability  analysis 
indicated  that  this  would  be  an  unstable  solution  (figure  9).  The  reason  for  this 
disagreement  in  the  results  is  investigated  by  adding  an  initial  angle  of  roll  to  the  numerical 
integration  analysis.  Adding  a  small  angle  of  initial  roll  (<p0  =  1  degree)  caused  the  vehicle 

to  steady  out  at  137  degrees  vice  159  degrees  as  shown  in  figure  34.  This  initial  roll  angle 
also  caused  a  steady  state  roll  angle  of  17  degrees  as  shown  in  figure  35.  In  turn,  this 
steadv  state  roll  angle  caused  the  steadv  state  yaw  velocity  (r)  shown  in  figure  36  (i.e. 
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motion  is  no  longer  restricted  to  the  vertical  plane).    Figure  37  shows  a  plot  of  z  versus 
and  indicates  that  the  vehicle  is  taking  a  helical  ascent  as  dicussed  by  Booth  [Ref  2:  304 
305].  Therefore,  the  numerical  integration  solution  resulting  in  a  steady  state  pitch  angle  c 
159  degrees  (figure  33)  is  unstable  in  the  horizontal  plane  as  predicted  by  the  linearize 
dynamic  response  analysis. 
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Figure  33.      Numerical   Integration   Solution   for  Angle  of  Pitch   (0)   when 
Center  of  Gravity  is  Aft  (xQB  =  -  1%) 
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Figure  34.     Numerical  Integration   Solution  for  Angle  of  Pitch   (0)   when 
xrD  =  -  \°7c  and  Initial  Roll  Angle  (<p, )  is  1   degree 
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Figure  35.       Numerical   Integration  Solution  for  Angle  of  Roll  (<p)   when 
XGB  =  "  J^c  and  ^mtia'  R°H  Angle  (<p0)  is  1  degree 
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Figure  36.       Numerical  Integration  Solution  for  Yaw  Velocity  (r)  when 
XGB  =  "  *^f  and  I™1'3'  Roll  Angle  (<p0)  is  1  degree 
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Figure  37.     Numerical  Integration  Solution  for  Heave  Velocity  (z)  when 
XGB  =  *  *^  anc*  ^mt'al  R°H  Angle  (<p0)  is  1  degree 
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Some  question  still  remains  with  regards  to  the  measure  of  stabilty  of  the  'inverted 

pendulum'  solutions  predicted  by  both  the  linearized  dynamic  response  analysis  and  the 

numerical  integration  analysis.   The  linearized  dynamic  response  analysis  predicts  a  stable 

solution  for  the  case  when  center  of  gravity   is   placed   aft   of  center   of  buoyancy 
(xGB  =  - 1  %)  and  the  dive  plane  angle  (&s)  is  -  7  degrees  (figure  10).    The  corresponding 

steady  state  value  for  pitch  angle  was  118  degrees  (figure  7).    A  random  persistent  roll 
disturbance  (<p ,)  was  added  to  the  numerical  integration  model  and  the  results  are  shown  in 

figure  38.    The  solid  line  indicates  the  results  when  a  small  disturbance  is  added  ($ 

centered  about  0.1  degrees),  the  dashed  line  indicates  the  results  when  a  large  disturbance 
is  added  (<j>,  centered  about  1.0  degrees).    As  expected,  the  large  disturbance  caused  the 

vehicle  to  roll  over  as  shown  by  the  resulting  angles  of  roll  ((p)  in  figure  39.   However,  the 

vehicle  continued  to  remain  stable  during  small  constant  random  disturbances  in  the 

inverted  position  as  shown  by  the  resulting  angles  of  roll  (<{>)   in  figure  40.   This  indicates 

that  indeed  these  'inverted  pendulum'  solutions  have  a  significant  measure  of  stability. 
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Figure  38.     Numerical  Integration   Solution  for  Angle  of  Pitch  (0)   when 


GB 


1%  and  a  Persistent  Roll  Disturbance  is  Added 
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Figure  39.     Numerical  Integration  Solution  for  Angle  of  Roll  (<p)  when 
x„B  =  -  1%   and  a  Large  Persistent  Roll  Disturbance  is  Added 
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Figure  40.     Numerical   Integration  Solution  for  Angle  of  Roll  (<p)   when 


GB 


Wc   and  a  Small  Persistent  Roll  Disturbance  is  Added 
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VI.    CONCLUSIONS  AND  RECOMMENDATIONS. 

•  The  steady  state  analysis  resulted  in  four  possible  solutions  provided  that  vehicle 
motion  was  restricted  to  the  vertical  plane.  Analyzing  the  dynamic  stability  using  the 
steady  state  results  as  nominal  points  generally  indicates  that  only  one  (if  any)  of  the 
four  possible  solutions  will  be  stable.  There  are  a  few  cases  where  two  solutions  are 
stable,  but  these  cases  (certain  non-zero  bow  plane  deflection  angles)  appear  to  be  the 
exception  and  not  the  rule. 

•  The  dynamic  stability  characteristics  of  submersibles  can  be  separated  with  respect  to 
vertical  plane  motions  (u,w,q,0)  and  horizontal  plane  motions  (v,p,r,(f)). 

•  It  is  possible  for  submersibles  to  be  dynamically  stable  with  respect  to  vertical  plane 
motion  in  the  inverted  (belly  up)  position  during  ascents  ('Inverted  Pendulum' 
stabilization). 

•  'Inverted  Pendulum'  stabilization  is  also  possible  in  the  horizontal  plane. 

•  Submersibles  are  able  to  maintain  this  inverted  orientation  (i.e.  ascend  belly  up 
without  rolling  over)  even  under  some  persistent  roll  excitation. 

•  As  a  recommendation,  the  dynamic  stability  analysis  should  be  expanded  to  include 
the  case  where  angle  of  roll  is  180  degrees,  and  the  case  where  angle  of  roll  is  non- 
zero (i.e.  <|)  neither  equals  zero  nor  180  degrees).  Analyzing  the  (p  =180  degrees  cases 
will  only  involve  changing  a  few  signs  with  regards  to  trigonometric  functions; 
however,  analyzing  the  non-zero  cases  will  require  significant  effort. 

•  Furthermore,  identifying  and  characterizing  different  stability  regions  over  a  range  of 
variations  of  the  system  parameters  should  be  the  matter  of  future  research. 
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APPENDIX  A 

SIX  DEGREE  OF  FREEDOM  EQUATIONS  OF  MOTION 

Source:  Smith,  Crane,  and  Summey  [Reference  1:11-16] 

1.  LONGITUDINAL  (SURGE)   EQUATION   OF  MOTION 

m  i  u  -  vr  +  wq  -  xq  (q2  +  r2)  +  yG  (pq  -  r)  +  zG  (pr  +  q):  = 

[  Xpp  p2  +  Xqq  q2  +  Xrr  r2  +  Xpr  pr 
+  Xu  u  +  XWq  wq  +  XVp  vp  +  Xvr  vr 
+  uq  (  Xq6s  6S  +  Xq6b  6b)  +  Xr5r  ur6r 
+    Xvv  v2  +  Xww  w2  +  Xv6r  uv6r  +  uw  (  XW6S  6S  +  Xwftb  6b) 

+  u2  (  X6S6S  6S2  +  X6b6b  6b2  +  Xar6r  6r2)  -  (W  -  B)  sin  9 

+  11-  Xprop 

2.  LATERAL  (SWAY)   EQUATION   OF  MOTION 

m  v  +  ur  -  wp  +  xq  (pq  +  r)  -  \q  (p2  +r2)  +  zq  (qr  -  p)  = 
Ypp  + Yr  r  +  Ypqpq+ Yqr  qr^ 
+    Yv  v  +  Ypup  +  Yr  ur  +  Yvq  vq  +  Ywp  wp  +  Ywr  wr 

+    Yv  uv  +  YVw  vw  +  Y&r  u-  6r 

'xnose  (    ,  .  \ 

CDy  h(x)  (v+xr-j2  +CDz  b(x)  (w-xq)2  j  ™  dx 

Atail 
+  (W-B)  cos  9  siiKJ) 


1} 


NORMAL  (HEAVE)  EQUATION  OF  MOTION 

m  [w  -  uq  +  vp  +  xG  (pr  -  q)  +  yG  (qr  +    p)  -  zG  (p2  +  q2)'  = 
[  Zqq  +  Zpp  p2  +  Zpr  pr  +  Zrr  r2] 
+  [Zw  w  +  Zquq  +  ZVp  vp  +  Zvr  vr  ] 
+  [Zw  uw  +  Zvv  v2  +  u2  (Zss  6S+Zfcb  6b)j 


'nose 


CDy  h(x)  (v+xr)2  +CDz  b(x)  (w-xq)2  ]  <^j  dx 


xtail 
+  (W-B)cos6  cos  (J) 

ROLL  EQUATION  OF  MOTION 

h  P  +  (Iz  -lyl)  qr  +  lxy  (Pr  "  q) "  Iyz  (q2  ■ r2) 

-  Ixz  (pq  +  r )  +  m  YG  (w  -  uq  +  vp)  -  zq  (v  +  ur  -  wp)  = 
JKp  p  +Kr  r   +  Kpq  pq  +  Kqr  qr 

+  ;  Kv  v    +  Kp  up  +  Kr  ur  +  KVq  vq  +  KWp  wp  +  Kwr  wr 
+  [  Kv  uv  +  Kvw  vwj  +  (yG  W  -  yg  B)  cos  6  cos  <p 
-  (zGW  -  zgB)  cos  0  sin  (j)  +  u^  Kprop 
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5.  PITCH  EQUATION  OF  MOTION 

Iy  q  +  (Ix  -  Iy)  pr  -  Ixy  (qr  +  p)  +  Iyz  (pq  -  r) 

+  Ixz  (p2  -  r2)  -  m  xq  (w  -  uq  +  vp)  -  zq  (u  -  vr  +  wp)  = 
^Mq  q   +  Mpp  p2  +  Mpr  pr  +  Mrr  r2] 
+  Mww   +  Mq  uq  +  MVp  vp  +  Mvr  vr  ] 

+  I  Mw  uw  +  Mvv  v2  +   u2  (M5S  6S  +  M5t>  6b),  ■ 

/"xnose  .  ,         . 

+  CDy  h(x)  (v+xr)2  +CDz  b(x)  (w-xq)2  I  Bl  x  dx 

ix,   •,  "  Uct<X> 

-  (xq  W  -  xg  B)  cos  6  cos  o  -  (zq  W  -  zg  B)  sin  6 

6.  YAW  EQUATION  OF  MOTION 

Iz  r  +  (Iy  -  Ix)  pq  -  Ixy  (p2  -  q2)  -  Iyz  (pr  +  q) 

+  Ixz  (qr  -  p)  -  m  xG  (v  +  ur  -  wp)  -  yG  (u  -  vr  +  wq)  = 
[Np  p  +Nr  r    +  Npq  pq  +  Nqr  qr 
+  Nvv   +  Np  up  +   Nr  ur  +  NVq  vq  +NWn  wp  +  Nwr  wr 

+    Nv  uv  +  Nvw  vw  +  N5r  u2  6r 

'xnose-  {  +  ■  \ 

CDy  h(x)  (v+xr)2  +  CDz  b(x)  (w-xq)2    [-^\  x  dx 

xtail  UdW 

+  (xq  W  -  xg  B|  cos  6  sin  <}>  +  (vq  W  -  yg  B)  sin  6  +  u2  Nprop 
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APPENDIX  B 

ROTATION  SEQUENCE  AND  EULER  ANGLE  RATES 

1 .  ROTATION  SEQUENCE  FOR  <j>,  0  AND  ^ 

Smith,  Crane,  and  Summey  [Reference  1:8]  descibe  the  transition  from  body  fixed  to 

inertia]  reference  frames  as  follows: 

Since  the  equations  of  motion  are  referred  to  an  axis  system  that  is  fixed  for  the 
SDV  (swimmer  delivery  vehicle),  and  thus  translates  and  rotates  with  it,  the 
orientation  and  position  of  the  moving  body  axis  system  relative  to  a  fixed  inertial 
reference  system  must  be  specified.  The  orientation  of  the  body  axis  system  with 
respect  to  the  inertial  reference  system  is  defined  by  the  standard  Euler  angles  V 
(yaw),  0  (pitch),  and  §  (roll).    The  rotation  sequence  from  the  inertial  reference 

system  to  the  body  system  is  H'  ,0  ,  and  <j>  as  shown  in  Figure  Bl  taken  from  Smith, 
Crane,  and  Summey  [Reference  1:18]. 

2.  EULER  ANGLE  RATES  FOR  <\>,Q  AND  V 

The  Euler  angle  rates  used  along  with  the  six  equations  of  motion  (Appendix  A)  in 
order  to  completely  determine  the  motion  of  the  submersible  were  specified  by  Smith, 
Crane,  and  Summey  [Reference  1:20]  to  be: 

<{>    =  p  +  q  sin  ty  tan  0  +  r  cos  <})  tan  0 

0    =  q  cos  (p  -  r  sin  (j) 

sin  0        cos  (b 

xp    =  q 31  +  1 

cos  0        cos  0 
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(l)  Vehicle-Centered 
Gravity-Directed  System 
para  1 le 1  to  inert  i  a  1 
axis  system. 


(2)  Horizontal  Flight 
Reference  System  derived 
from  XqYqZq  by  rotation  a- 
bout  Z  through  yaw  angle  i> . 


Y*  Y* 


(3)  Rol 1-Resolved  F)  ight 
Reference  System  derived 
from  X"Y"Z"  by  rotation 
about  Y"  through  pitch 
ang  le  9 . 


(X1) 

(h)  Vehicle  Body  Axis  Ref- 
erence System  derived  from 
X'Y'Z'  by  rotation  about  X 
through  rol 1  angle  $ . 


Figure  Bl.     Unit  Sphere  Development  of  Euler  Angles 
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STEADY  STATE  COMPUTER  PROGRAM 

PROGRAM       STEADY 

STEADY    STATE    SOLUTIONS    IN    THE    VERTICAL    PLANE 

DIVE    PLANE    VARIATION 

REAL  L,MASS,IX,IY,IZ,IXZ,IYZ,IXY,LAMBDA 

REAL  KPDOT , KRDOT , KPQ , KQR , KVDOT , KP , KR , KVQ , KWP , KWR , KV 

REAL  KVW,KPN,KDB 

REAL  MQDOT , MPP , MPR , MRR , MWDOT , MQ , MVP , MVR , MW , MW , MDS 

REAL  MDB,NDRB 

REAL  NPDOT , NRDOT , NPQ , NQR , NVDOT , NP , NR , NVQ , NWP , NWR , NV 

REAL  NVW,NDRS 

DIMENSION  X(9) ,BR(9) ,HH(9) ,VEC1(9) 

GEOMETRIC  PROPERTIES  AND  HYDRODYNAMIC  COEFFICIENTS 


PI 
WEIGHT= 

L 

RHO 

G 

CDO 

MASS 

CDZ 

xww 

XWDS 

XWDB 

XDSDS 

XDBDB 

CDO 

ZW 

ZDS 

ZDB 

MW 

MDS 

MDB 


A . 0*ATAN(1 .0 ) 
12000.0 

17.425 

1.94 
32.2 
0.0057 
=WEIGHT/G 
=0.5*0.5*RHO 


=  1 
=  4 
=  9 
=  -1 
=  -8 

=  -3 
=  -2 
=  -2 

=  9 
=  -1 
=  1 


710E-01*0 
600E-02*0 
660E-03*0 
160E-02*0 
070E-03*0 
CD0*0 
020E-01*0 
270E-02*0 
270E-02*0 
860E-02*0 
113E-02*0 
113E-02*0 


5*RHO* 
5*RHO* 
5*RHO* 
5*RHO* 
5*RHO* 
5*RHO* 
5*RHO* 
5*RHO* 
5*RHO* 
5*RHO* 
5*RHO* 
5*RHO* 


L**2 
L**2 
L**2 
L**2 
L**2 
L**2 
L**2 
L**2 
L**2 
L**3 
L**3 
L**3 


OPEN (21 ,NAME='ST1 .RES' , STATUS= ' NEW ' 
OPEN( 2  2,NAME='ST2.RES' , STATUS= ' NEW ' 
OPEN (23, NAME= ' ST3 . RES ' , STATUS= ' NEW' 
OPEN( 2  4 ,NAME='ST4 .RES' , STATUS= ' NEW ' 
OPEN ( 31 ,NAME='COEF.DAT' , STATUS =' NEW 


78 


105.9/12 

0 

-99.3/12 

0 

-87.3/12 

0 

-66.3/12 

0 

72.7/12 

0 

83.2/12 

0 

91.2/12 

0 

99.2/12 

0 

103.2/12 

0 

0.00/12 

0 

8.24/12 

0 

19.76/12 

0 

29.36/12 

0 

31.85/12. 

0 

27.84/12 

0 

21.44/12. 

0 

12.00/12. 

0 

0.00/12 

0 

C      DEFINE  THE  LENGTH  X,  BREADTH  BR,  AND  HEIGHT  HH  TERMS 
C 

X(l)   = 

X(2)   = 

X(3)   = 

X(4)   = 

X(5)   = 

X(6)   = 

X(7)   = 

X(8)   = 

X(9)   = 
C 

BR(1)  = 

BR(2)  = 

BR(3)  = 

BR(4)  = 

BR(5)  = 

BR(6)  = 

BR(7)  = 

BR(8)  = 

BR(9)  = 
C 

C      COMPUTE  AREA  AND  CENTROID 
C 

CALL  TRAP( 9,BR,X,AREA) 

DO  9  1=1,9 

VEC1 ( I )=X( I ) *BR( I ) 
9  CONTINUE 

CALL  TRAP(9,VEC1,X,XAA) 

XA=XAA/AREA 
C 

WRITE  (*,1002) 

READ   (*,*)      DSMIND,DSMAXD, IDS 

DSMIN=DSMIND*PI/18  0 

DSMAX=DSMAXD*PI/18  0 

WRITE  (*,1001) 

READ   (*,*)      RATIO 

WRITE  (*,1003) 

READ   (*,*)      DELB 

DELB=DELB*WEIGHT/100.0 

WRITE  (*,1004) 

READ   (*,*)     XGB 

XGB=XGB*L/100. 0 

WRITE  (*,1005) 

READ   ( * , * )     ZGB 

WRITE  (*,1006) 

READ   ( * , * )     XB 

XB=XB*L/100 . 0 

WRITE  (*,1007) 

READ   (*,*)      ZB 

WRITE  (31,*)   RATIO,  DELB,  XGB,  ZGB,  XB ,  ZB 
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DO  1  1=1, IDS 

DS=DSMIN+(DSMAX-DSMIN)*(I-1)/(IDS-1) 
IF  (DELB.EQ.0.0)  DELB=0 . 000001 
IF  (ZGB.EQ.0.0)   ZGB  =0.000001 
DB=RATIO*DS 

PX     =XGB*WEIGHT-XB*DELB 

PZ     =ZGB*WEIGHT-ZB*DELB 

DEN    =CDZ*AREA*(PX+XA*DELB) 

LAMBDA=MW*DELB-PX*ZW+PZ*(XWDS+RATIO*XWDB)*DS 

ALPHA  — PX*(ZDS  +  RATIO*ZDB)*DS-PZ*CD0  +  PZ*'(XDSDS  + 
&  RATIO*RATIO*XDBDB)*DS**2+DELB*(MDS+RATIO 

&  *MDB)*DS 

BETA   =PZ*XWW 

LAMBDA=LAMBDA/DEN 

ALPHA  =ALPHA  /DEN 


BETA   =BETA   /DEN 

c 

c 

A   =  1.0+BETA 

B   =  LAMBDA 

C   =  ALPHA 

DET=  B**2-4.0*A*C 

IF  (DET.LT.0.0)  GO 

TO  2 

WP=(-B+SQRT(DET) )/( 2.0*A) 
YY=-XWW*WP**2-(XWDS*DS+XWDB*DB)*WP 

-(XDSDS*DS**2+XDBDB*DB**2)+CD0 
IF  (WP.GE.0.0)  XX=ZW*WP+ZDS*DS+ZDB*DB-CDZ*AREA 

*WP*ABS(WP) 
IF  (WP.LT.0.0)  XX=ZW*WP+ZDS*DS+ZDB*DB+CDZ*AREA 

*WP*ABS(WP) 
THETA=ATAN2 ( YY , XX ) 
USQ=DELB*SIN(THETA)/YY 
THETA=THETA*18  0/PI 
DSD=DS*180/PI 
IF  (USQ.LT.0.0)  GO  TO  3 
IF  (WP.GE.0.0)  U=  SQRT(USQ) 
IF  (WP.LT.0.0)  U=-SQRT(USQ) 
W=WP*U 
WRITE  (21,*)  DSD,THETA,U,W,WP 

WP=(-B-SQRT(DET) )/(2.0*A) 
YY=-XWW*WP**2-(XWDS*DS+XWDB*DB) *WP 

-(XDSDS*DS**2+XDBDB*DB**2 )+CD0 
IF  (WP.GE.0.0)  XX=ZW*WP+ZDS*DS+ZDB*DB-CDZ*AREA 

*WP*ABS(WP) 
IF  (WP.LT.0.0)  XX=ZW*WP+ZDS*DS+ZDB*DB+CDZ*AREA 

*WP*ABS(WP) 
THETA=ATAN2 ( YY , XX ) 
USQ=DELB*SIN(THETA)/YY 
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DSD=DS*180/PI 

THETA=THETA*18  0/PI 

IF  (USQ.LT.O.O)  GO  TO  2 

IF  (WP.LT.0.0)  U=-SQRT(USQ) 

IF  (WP.GE.0.0)  U-  SQRT(USQ) 

W=WP*U 

WRITE  (22,*)  DSD,THETA,U,W,WP 


2    A   =  -1.0+BETA 

DET=  B**2-4.0*A*C 
IF  (DET.LT. 0.0)  GO  TO  1 
WP=(-B+SQRT(DET) )/(2.0*A) 
YY=-XWW*WP**2-(XWDS*DS+XWDB*DB)*WP 
&  -(XDSDS*DS**2+XDBDB*DB**2 )+CD0 

IF  (WP.LT.0.0)  XX=ZW*WP+ZDS*DS+ZDB*DB-CDZ*AREA 

*WP*ABS(WP) 
IF  (WP.GE.0.0)  XX=ZW*WP+ZDS*DS+ZDB*DB+CDZ*AREA 

*WP*ABS(WP) 
THETA=ATAN2 ( YY , XX ) 
USQ=DELB*SIN(THETA)/YY 
DSD=DS*180/PI 
THETA=THETA*18  0/PI 
IF  (USQ.LT. 0.0 )  GO  TO  4 
IF  (WP.GE.0.0)  U--SQRT(USQ) 
IF  (WP.LT.0.0)  U«  SQRT(USQ) 
W=WP*U 

WRITE  (23,*)  DSD,THETA,U,W,WP 
4    WP=(-B-SQRT(DET) )/(2.0*A) 

YY=-XWW*WP**2-(XWDS*DS+XWDB*DB ) *WP 
&  -(XDSDS*DS**2+XDBDB*DB**2 )+CD0 

IF  (WP.LT.0.0)  XX=ZW*WP+ZDS*DS+ZDB*DB-CDZ*AREA 

*WP*ABS(WP) 
IF  (WP.GE.0.0)  XX=ZW*WP+ZDS*DS+ZDB*DB+CDZ*AREA 

*WP*ABS(WP) 
THETA=ATAN2 ( YY , XX ) 
USQ=DELB*SIN(THETA)/YY 
DSD=DS*180/PI 
THETA=THETA*18  0/PI 
IF  (USQ.LT.O.O)  GO  TO  1 
IF  (WP.GE.0.0)  U—  SQRT(USQ) 
IF  (WP.LT.0.0)  U-  SQRT(USQ) 
W=WP*U 
WRITE  (24,*)  DSD,THETA,U,W,WP 

1  CONTINUE 

STOP 

1001  FORMAT  ( '  ENTER  BOW  PLANE  TO  DIVE  PLANE  RATIO' ) 

1002  FORMAT  ('  ENTER  MIN,  MAX,  AND  INCREMENTS  IN 
&  DS  (degrees)') 

1003  FORMAT  ('  ENTER  DELB  (  %W )  '  ) 

1004  FORMAT  ('  ENTER  XGB  (%L)') 


SI 


1005 

FORMAT  ('  ENTER  ZGB  (feet)') 

1006 

FORMAT  ('  ENTER  XB  ( %L ) ' ) 

1007 

FORMAT  ('  ENTER  ZB  (feet)') 

c 

END 

c 
c 

SUBROUTINE  TRAP ( N , A , B , OUT ) 

NUMERICAL  INTEGRATION  ROUTINE  USING 

c 
c 

THE  TRAPEZOIDAL  RULE 

DIMENSION  A(l) ,B(1) 

Nl=N-l 

OUT=0 . 0 

DO  1  1=1, Nl 

OUT1=0.5*(A(I )+A(I+l) )*(B(I+1)-B(I ) 

OUT  =OUT+OUTl 

1 

CONTINUE 

RETURN 

END 
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APPENDIX  D 


LINEARIZED  DYNAMIC  STABILITY  COMPUTER  PROGRAM 

C      PROGRAM  LINEARIZED  DYNAMIC  STABILITY 

C       10        20        30         40         50 

C2 34 56 78901234567890123456789012345678901234567890123456 

C 

IMPLICIT  DOUBLE  PRECISION  (A-H,0-Z) 

DOUBLE  PRECISION  L , MASS , IA , IX , I Y , IZ 

DOUBLE  PRECISION  KPDOT , KRDOT , KPQ , KQR , KVDOT , KP , KR , 
&   KVQ , KWP , KWR , KV , 

&   KW , KPN , KDB , MQDOT , MPP , MPR , MRR , MWDOT , MQ , MVP , MVR , 
&   MW,MW,MDS,MDB, 

&   NDRE , NPDOT , NRDOT , NPQ , NQR , NVDOT , NP , NR , NVQ , NWP , NWR , 
&   NV,NVW,NDRS 
C 

DIMENSION  Al(4,4),Bl(4,4) ,BETAl( 4 ) ,ALFR1 ( 4 ) ,ALFI1 ( 4 ) 

DIMENSION  BBl(4,4)/BB2(4,4),ZZZl(4,4),ZZZ2(4,4) 

DIMENSION  A2(4,4),B2(4,4) , BETA2(4 ) ,ALFR2( 4 ) ,ALFI2( 4 ) 

DIMENSION  WR1 ( 4 ) ,WR2 ( 4 ) ,WIl ( 4 ) ,WI2 ( 4 ) 

DIMENSION  X( 9 ) ,BR(  9 ) ,VEC1 ( 9 ) ,VEC2 ( 9 ) 
C 

C      GEOMETRIC  PROPERTIES 
C 

PI     =    4 .D0*DATAN(1 .DO ) 

WEIGHT-    12000.0 

IX     =  1760.0 

IY     =  9450.0 

IZ     =10700.0 

L      =    17.425 

RHO    =    1.94 

G      =    32.2 

CD0    =    0.0057 

MASS   =    WEIGHT/G 

CDZ    =    0.5*0.5*RHO 

CD0    =    CD0*0 . 5*RHO*L**2 
C 

C      SURGE  HYDRODYNAMIC  COEFFICIENTS 
C 


XPP    =  7 

.030E- 

-03*0 

5*RHO'L**4 

XQQ    =-1 

.470E- 

-02*0 

5*RHO*L**4 

XRR    =  4 

.010E- 

-03*0 

5*RHO*L**4 

XPR    =  7 

.640E- 

-04*0 

5*RHO*L**4 

XUDOT  =-7 

.  580E- 

-03*0 

5*RHO*L**3 

XWQ    =-1 

.920E- 

-01*0 

5*RHO*L**3 
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c 
c 
c 


c 
c 

c 


XVP 

=  - 

-3 

.240E- 

-03*0. 

5*RHO*L**3 

XVE 

= 

1 

.890E- 

-02*0. 

,  5*RHO*L**3 

XQDS 

= 

2 

.610E- 

-02*0. 

5*RHO*L**3 

XQDB 

=  - 

-2 

.600E- 

-03*0. 

5*RHO*L**3 

XRDR 

=  - 

-8 

.180E- 

-04*0. 

,  5*RHO*L**3 

XW 

= 

5 

.290E- 

-02*0. 

,  5*RHO*L**2 

XWW 

= 

1 

.710E- 

-01*0. 

,  5*RHO*L**2 

XVDR 

= 

1 

.730E- 

-03*0. 

,5*RHO*L**2 

XWDS 

= 

4 

.600E- 

-02*0. 

,  5*RHO*L**2 

XWDB 

= 

9 

.660E- 

-03*0. 

,  5*RHO*L**2 

XDSDS 

=  - 

-1 

.160E- 

-02*0. 

,  5*RHO*L**2 

XDBDB 

=  - 

-8 

.070E- 

-03*0. 

,  5*RHO*L**2 

XDRDR 

=  - 

-1 

.010E- 

-02*0. 

,  5*RHO*L**2 

XRES 

= 

CDO*0. 

,5*RHO*L**2 

SWAY  HYDRODYNAMIC  COEFFICIENTS 


YPDOT 

= 

1, 

.270E- 

-04*0. 

, 5*RHO*L**4 

YRDOT 

= 

1, 

.240E- 

-03*0. 

, 5*RHO*L**4 

YPQ 

= 

4, 

.125E- 

-03*0, 

,5*RHO*L**4 

YQR 

=  • 

-6, 

.510E- 

-03*0. 

, 5*RHO*L**4 

YVDOT 

=  ■ 

-5, 

.550E- 

-02*0, 

. 5*RHO*L**3 

YP 

= 

3, 

.055E- 

-03*0, 

, 5*RHO*L**3 

YR 

= 

2, 

.970E- 

-02*0, 

,5*RHO*L**3 

YVQ 

= 

2. 

.360E- 

-02*0. 

,5*RHO*L**3 

YWP 

= 

2, 

.350E- 

-01*0, 

,5*RHO*L**3 

YWR 

=  - 

-1, 

.880E- 

-02*0, 

, 5*RHO*L**3 

YV 

=  - 

-9, 

.310E- 

-02*0, 

, 5*RHO*L**2 

YWJ 

= 

6, 

.840E- 

-02*0, 

, 5*RHO*L**2 

YDRS 

=  +  2, 

.270E- 

-02*0, 

,5*RHO*L**2 

YDRB 

=  +  2, 

.270E- 

-02*0, 

, 5*RHO*L**2 

HEAVE  HYDRODYNAMIC  COEFFICIENTS 


ZQDOT 

=  ■ 

-6, 

.810E- 

-03*0, 

, 5*RHO*L**4 

ZPP 

= 

1, 

.270E- 

-04*0, 

.5*RHO*L**4 

ZPR 

= 

6, 

.670E- 

-03*0, 

, 5*RHO*L**4 

ZRR 

=  - 

-7, 

.350E- 

-03*0, 

,5*RHO*L**4 

ZWDOT 

=  ■ 

-2, 

.430E- 

-01*0, 

.  5*RHO*L**3 

ZQ 

=  - 

-1. 

.350E- 

-01*0, 

,5*RHO*L**3 

ZVP 

=  • 

-4 

.810E- 

-02*0, 

, 5*RHO*L**3 

ZVR 

= 

4, 

.550E- 

-02*0. 

,  5*RHO*L**3 

ZW 

=  - 

-3, 

.020E- 

-01*0, 

,  5*RHO*L**2 

ZVV 

=  - 

-6, 

.840E- 

-02*0, 

,  5*RHO*L**2 

ZDS 

=  - 

-2, 

.270E- 

-02*0, 

.  5*RHO*L**2 

ZDB 

=  - 

-2, 

.270E- 

-02*0, 

,  5*RHO*L*^2 

ROLL  HYDRODYNAMIC  COEFFICIENTS 


KPDOT  =-1 
KRDOT  =-3 
KPQ    =-6 


010E-03*0 
370E-05*0 
930E-05*0 


5*RHO*L**5 
5*RHO*L**5 
5*RHO*L**5 
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c 
c 
c 


c 

c 

c 


c 
c 

c 


KQR 

KVDOT 

KP 

KR 

KVQ 

KWP 

KWR 

KV 

KVW 


680E-02*0 
270E-04*0 
100E-02*0 
410E-04*0 
115E-03*0 
270E-04*0 
390E-02*0 
055E-03*0 
870E-01*0 


5*RHO* 
5*RHO* 
5*RHO* 
5*RHO* 
5*RHO* 
5*RHO* 
5*RHO* 
5*RHO* 
5*RHO* 


L**5 
L**4 
L**4 
L**4 
L**4 
L**4 
L**4 
L**3 
L**3 


PITCH  HYDRODYNAMIC  COEFFICIENT 


MQDOT 

MPP 

MPR 

MRR 

MWDOT 

MQ 

KVP 

MVR 

MW 

KW 

MDS 

MDB 


680E- 
260E- 
040E- 
860E- 
810E- 
860E- 
180E- 
730E- 
860E- 
510E- 
113E- 
113E- 


02*0 
05*0 
03*0 
03*0 

02*0 
02*0 
03*0 
02*0 
02*0 
02*0 
02*0 
02*0 


5*RHO* 
5*RHO* 
5*RHO* 
5*RHO* 
5*RHO* 
5*RHO* 
5*RHO* 
5*RHO* 
5*RHO* 
5*RHO* 
5*RHO* 
5*RHO* 


L**5 
L**5 
L**5 
L**5 
L**4 
L**4 
L**4 
L**4 
L**3 
L**3 
L**3 
L**3 


YAW  HYDRODYNAMIC  COEFFICIENTS 


NPDOT 

=  -3 

,370E-05*0. 5*RHO*L**5 

NRDOT 

=  -3, 

.  400 E- 03*0. 5*RHO*L**5 

NPQ 

=  -2, 

.  110E-02*0. 5*RHO*L**5 

NQR 

=  2, 

.750 E- 03*0. 5*RHO*L**5 

NVDOT 

=  1. 

.  240 E- 03*0. 5*RHO*L**4 

NP 

=  -8, 

.  4  05E-04*0 . 5*RHO*L**4 

NR 

—  1. 

.  640 E- 02*0. 5*RHO*L**4 

NVQ 

=  -9, 

.  990 E- 03*0. 5*RHO*L**4 

NWP 

=  -1  , 

.  750 E- 02*0. 5*RHO*L**4 

NWR 

=  7  , 

.  350 E- 03*0. 5*RHO*L**4 

NV 

=  -7. 

,  420 E- 03*0. 5*RHO*L**3 

NVW 

=  -2, 

.670 E- 02*0. 5*RHO*L**3 

NDRS 

—  1, 

.  113E-02*0.5*RHO*L**3 

NDRB 

■+1, 

.113 E- 02*0. 5*RHO*L**3 

DEFINE  THE  LENGTH  X  AND  BREADTH 

X(l) 

= 

-105.9/12.0 

X(2) 

= 

-99.3/12.0 

X(3) 

= 

-87.3/12.0 

X(4) 

= 

-66.3/12.0 

X(5) 

= 

72.7/12.0 

X(6) 

= 

83.2/12.0 

X(7) 

= 

91.2/12.0 

X(8) 

= 

99.2/12.0 

BR    TERMS 


X(9 


103.2/12.0 
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BR(1) 

BR(2) 

BR(3) 

BR(4) 

BR(5) 

BR(6) 

BR(7) 

BR(8) 

BR(9) 

0 

.00/12, 

,0 

8 

.24/12. 

,  0 

19, 

.76/12. 

,0 

29, 

.36/12. 

,0 

31, 

.85/12. 

,  0 

27, 

.84/12. 

,0 

21, 

.44/12. 

.0 

12, 

.00/12. 

,0 

0, 

.00/12. 

,0 

c 

C      COMPUTE  AREA,  CENTROID,  AND  MOMENT  OF  INERTIA 
C 

CALL  TRAP( 9 , BR, X, AREA) 

DO  9  1=1,9 

VEC1 ( I )=X(I )*BR( I ) 
VEC2( I )=X(I )*VEC1(I ) 
9  CONTINUE 

CALL  TRAP( 9 ,VEC1 , X , XAA ) 

XA=XAA/AREA 

CALL  TRAP( 9,VEC2,X, IA) 
C 

WRITE  (*,1001) 

READ   (*,*)  IRES 

OPEN( 31 , NAME='COEF.DAT' , STATUS= ' OLD ' ) 

READ(31,*)   RATIO,  DELB,  XGB ,  ZGB,  XB ,  ZB 

BUOY=  WEIGHT  +  DELB 

XG=XB+XGB 

ZG=ZB+ZGB 
C 

C      MASS  MATRIX  COEFFICIENTS 
C 

Bl(l,l)=  MASS  -  XUDOT 

Bl(l, 2)=  0.0 

Bl( 1 , 3 )=  MASS*ZG 

Bl(l, 4)=  0.0 
C 

Bl(2,l)=  0.0 

Bl(2,2)=  MASS  -  ZWDOT 

Bl(2, 3  )=-(ZQDOT+MASS*XG) 

Bl(2,4)=  0.0 
C 

Bl( 3,1  )=  MASS*ZG 

Bl ( 3 , 2 ) =- ( MWDOT+MASS*XG ) 

Bl(3,3)=  IY-MQDOT 

Bl( 3,4  )=  0.0 
C 

Bl(4,l)=  0.0 

Bl(4,2)=  0.0 

Bl(4,3)=  0.0 

Bl(4,4  )=  1.0 
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B2( 1 , 1 )=  IX-KPDOT 

B2(l,2)=  0.0 

B2( 1 , 3 )=-( KVDOT+MASS*ZG) 

B2( 1 , 4  )=-KRDOT 
C 

B2(2,l  )=  0.0 

B2(2,2)=  1.0 

B2(2,3)=  0.0 

B2(2,4  )  =  0.0 
C 

B2( 3,1  )=-( YPDOT+MASS*ZG) 

B2(3,2)=  0.0 

B2(3,3)=  MASS-YVDOT 

B2(3,4)=  MASS*XG-YRDOT 
C 

B2 ( 4 , 1  )=-NPDOT 

B2(4,2  )  =  0.0 

B2(4,3)=  MASS*XG-NVDOT 

B2(4,4)=  IZ-NRDOT 
C 

OPEN( 41 , NAME='DEOS .RES' , STATUS= ' NEW ' 

OPEN(4  2,NAME=fDEOSl .RES' , STATUS= ' NEW 

OPEN( 43,NAME='DEOS2.RES' , STATUS* 'NEW ) 
C 

IF  ( IRES.EQ.l )  GO  TO  1 

IF  ( IRES.EQ.2 )  GO  TO  2 

IF  (IRES.EQ.3)  GO  TO  3 

IF  ( IRES.EQ. 4 )  GO  TO  4 

1  OPEN  ( 21 ,NAME='STl .RES' , STATUS= ' OLD ' ) 

11  READ  ( 21 , * ,END=100 )  DSD, THET0 ,U0 ,W0 , WP 
GO  TO  5 

2  OPEN  ( 22,NAME=' ST2 .RES' ,STATUS= 'OLD' ) 

12  READ  ( 22 , * ,END=100 )  DSD , THET0 , U0 ,W0 , WP 
GO  TO  5 

3  OPEN  ( 2  3 ,NAME='ST3.RES' , STATUS- ' OLD ' ) 

13  READ  ( 23, * ,END=100 )  DSD , THET0 , U0 ,W0 , WP 
GO  TO  5 

4  OPEN  ( 24 ,NAME='ST4 .RES' ,STATUS= 'OLD' ) 

14  READ  ( 24 ,* , END=100 )  DSD , THET0 , U0 , W0 , WP 
GO  TO  5 

C 

5  THETA0=THET0*PI/180. 0 
DS  =  DSD*PI/180.0 

DB  =  DS  *  RATIO 
C 

C      DAMPING  MATRIX  COEFFICIENTS 
C 

Al ( 1 , 1  )=-2 . 0*U0*CD0+W0* (XWDS*DS+XWDB*DB ) 
&  +2 . 0*U0* (XDSDS*DS**2+XDBDB*DB**2 ) 

Al  (1,2  )  =  2 . 0*XWW*WO  +  U0* ( XWDS*DS  +  XWDB*DB 
Al  (  1  ,  3  )=  (XWQ-MASS)*W0-(XQDS^D?^XQDB*DB)  *U0 
Al (1 , 4 )=-(WEIGHT-BUOY)-DCOS(THETA0 ) 
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c 
c 
c 


Al 

.1  ) 

Al 

2 

,2) 

Al 

2 

r3) 

Al 

"5 

.4  ) 

Al 

3 

,1  ) 

Al 

.3 

-2) 

Al 

3 

,3) 

Al 

3 

r  4  ) 

Al 

4 

-1) 

Al 

4 

-2) 

Al 

,4 

-3) 

Al 

4 

-4) 

A2 

1 

1) 

A2 

1 

-2) 

A2 

1 

,3) 

A2 

1 

-4) 

A2 

2 

1) 

A2 

2 

2) 

A2 

2 

3) 

A2l 

2 

4  ) 

A2( 

3 

1  ) 

A2! 

3 

2) 

A2I 

3 

3) 

A2I 

3 

4  ) 

A2I 

4 

1  ) 

A2I 

4 

2) 

A2l 

4 

3) 

A2l 

4 

4  ) 

ZW*W0+2 . 0*U0* (ZDS*DS+ZDB*DB ) 
ZW*U0-2. 0*CDZ*AREA*DABS(WO ) 
( ZQ+MASS ) *U0+2 . 0*CDZ*AREA*XA*DABS (WO ) 
•  (WEIGHT-BUOY ) *DS IN ( THETAO ) 

MW*W0+2 .0*UO* (MDS*DS+MDB*DB) 
MW*U0+2 . 0*CDZ*AREA*XA*DABS (WO ) 
(MQ-MASS*XG) *U0-MASS*ZG*W0 

-2.0*CDZ*IA*DABS(W0) 
( XG*WEI GHT-XB*BUOY ) *DSIN ( THETAO ) - 
( ZG*WEIGHT-ZB*BUOY) *DCOS(THETAO ) 

0.0 
0.0 
1.0 
0.0 

KP*UO+(KWP-MASS*ZG)*WO 
■  ( ZG*WEIGHT-ZB*BUOY ) *DCOS ( THETAO ) 
KV*UO+KVW*WO 
( KR+MASS*ZG ) *UO+KWR*WO 

1.0 
0.0 
0.0 
DTAN( THETAO) 

YP*UO+( YWP+WASS ) *W0 
(WEIGHT-BUOY ) *DCOS ( THETAO ) 
YV*UO+YVW*WO-CDZ*AREA*DABS (WO ) 
YWR*WO+ ( YR-MASS ) *UO-CDZ*AREA*XA 
*DABS(WO ) 

MASS*XG*WO+NP*UO+NWP*WO 

( XG*WEIGHT-XB*BUOY ) *DCOS ( THETAO ) 

NV*UO+NVW*WO-CDZ*AREA*XA*DABS(WO) 
(NR-MASS*XG)*UO+NWR*WO-CDZ*IA*DABS(WO ) 


RESTORE  B-MATRIX 


DO  71  1=1,4 

DO  72  J=l,4 

BB1( I, J)-B1(I , J 

72 

CONTINUE 

71 

CONTINUE 

DO  81  1=1,4 

DO  82  J=l,4 

BB2( I , J)=B2( I , J) 

CONTINUE 
CONTINUE 
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c 
c 


CALL  RGG( 4 , 4 ,Al ,BBl ,ALFRl ,ALFI1 ,BETAl , 0 ,  ZZZl , IER) 
CALL  DEGSTB(DE0S1 ,ALFR1 ,ALFI1 ,BETAl , FREQ1 ,WRl ,WI1 

CALL  RGG( 4 , 4 ,A2 ,BB2 ,ALFR2 ,ALFI2 , BETA2 , 0 ,ZZZ2 , IER) 
CALL  DEGSTB(DEOS2,ALFR2 ,ALFI2,BETA2 , FREQ2 ,WR2,WI2 

IF  (DE0S1 .GE.DE0S2 )  DEOS=DEOSl 
IF  (DEOS1 .LT.DEOS2)  DEOS-DEOS2 

WRITE  (41,2001)  DSD,THET0,U0,W0,WP,DEOS, 
&  DEOSl,DEOS2 

IF  (DEOS.LT. 0 .DO ) 
&    WRITE  (42,2001)  DSD,THET0,U0,W0,WP,DEOS, 
&  DEOSl,DEOS2 

IF  (DEOS1 .LT.0.D0) 
&    WRITE  (43,2001)  DSD , THET0 , UO , WO ,WP , DEOS , 
&  DEOSl,DEOS2 

IF  ( IRES.EQ.l )  GO  TO  11 

IF  ( IRES.EQ.2  )  GO  TO  12 

IF  ( IRES.EQ. 3)  GO  TO  13 

IF  ( IRES.EQ. 4  )  GO  TO  14 

i 

100  STOP 

1001  FORMAT  ('  ENTER  THE  RESPONSE  DATA  FILE  DESIRED 

&  (1,2,3,  OR  4)  ' ) 

2  0  01  FORMAT  (8E15.5) 

2002  FORMAT  (F10.3) 
END 

SUBROUTINE  DEGSTB ( DEOS , ALFR , ALFI , BETA , OMEGA , WR , WI ) 
IMPLICIT  DOUBLE  PRECISION  (A-H,0-Z) 
DIMENSION  ALFR( 4 ) ,ALFI ( 4 ) ,BETA( 4 ) ,WR( 4 ) ,WI ( 4 ) 
DO  1  1=1 , 4 

WR( I )=ALFR( I )/BETA( I ) 

WI ( I )=ALFI ( I )/BETA( I ) 

1  CONTINUE 
DEOS=-l .0E+10 
DO  2  1=1,4 

IF  (WR( I ) .LT.DEOS )  GO  TO  2 
DEOS=WR( I ) 
IJ  =  I 

2  CONTINUE 
OMEGA=WI ( I  J ) 
OMEGA=DABS ( OMEGA ) 
RETURN 

END 

SUBROUTINE  TRAP ( N , A , B , OUT ) 
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c 

C      NUMERICAL  INTEGRATION  ROUTINE  USING 

C  THE  TRAPEZOIDAL  RULE 

C 

IMPLICIT  DOUBLE  PRECISION  (A-H,0-Z) 

DIMENSION  A( 1 ) ,B(1 ) 

Nl-N-1 

OUT=0 . 0 

DO  1  1=1, Nl 

OUT1=0.5*(A(I )+A(I+l) )*(B(I+1 )-B(I ) 
OUT  =OUT+OUTl 
1  CONTINUE 
RETURN 
END 
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